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Abstract 

The  Air  Force  is  currently  developing  new  prodiicts  that  incorporate  a  variety 
of  chemicals  which  may  come  in  contact  with  prodiict  users.  To  define  which 
chemicals  are  dangerous  to  the  user,  toxicity  studies  have  been  performed.  However, 
analysis  of  toxicity  ultimately  requires  models  of  the  exposed  cellular  systems.  This 
thesis  provides  an  introduction  of  how  to  model  and  analyze  small  and  large  cellular 
systems.  Understanding  the  underlying  behavior  of  small  models  and  their  relation 
to  large  systems  will  lead  to  a  better  understanding  of  how  the  Air  Force  should 
construct  intracellular  models  to  assist  in  future  toxicology  studies. 

Developing  analysis  techniques  to  include  steady  state  analysis  through  lin¬ 
earization,  and  then  considering  small  reaction  systems,  such  as  the  Brusselator 
and  Schnackenberg  models,  led  to  a  basic  understanding  of  model  behavior.  This 
knowledge  was  applied  to  create  new  models  in  an  effort  to  begin  a  transition  from 
previously  created  models  to  the  construction  of  models  unique  to  the  Air  Force. 

Sensitivity  analyses  performed  on  existing  systems  furthered  research  efforts 
by  developing  knowledge  of  how  systems  behave  under  various  initial  conditions  and 
perturbations  of  uncertain  constant  parameters.  Analysis  displayed  great  sensitivity 
within  some  models.  This  analysis  was  applied  to  a  new  model  to  look  for  interesting 
behavior  such  as  oscillatory  convergence.  The  new  model  was  then  incorporated 
into  a  larger  model  to  determine  how  its  behavior  changed  with  respect  to  changes 
in  the  larger  model. 

This  knowledge  of  how  small  systems  behave  in  relation  to  larger  systems 
should  help  the  Air  Force  to  develop  and  analyze  intracellular  toxicology  models. 


Cell  Modeling 


I.  Introduction 

1.1  Overview 

“The  development  of  new  products  within  the  Air  Force  always  involves  the  use 
of  a  variety  of  chemicals  that  may  be  needed  in  the  prodiiction,  use,  or  maintenance 
of  the  product.  Often,  these  chemicals  are  either  novel  or  are  being  used  in  ways 
unique  to  the  Air  Force.  In  order  to  facilitate  the  development  of  these  prodiicts, 
while  maintaining  the  safety  of  Air  Force  personnel,  rapid  and  accurate  methods  for 
determining  the  toxicity  of  these  chemicals  are  needed.  Due  to  the  complexity  of 
human  biology,  standard  toxicological  assays  are  time  consuming  and,  due  to  the 
lack  of  adequate  understanding  of  toxicological  mechanisms,  may  be  inaccurate  at 
predicting  human  risk.  Because  of  these  limitations,  the  Air  Force  has  committed 
to  the  development  of  new  toxicity  assays  that  will  generate  large  amounts  of  data 
on  gene  expression  and  protein  concentration... The  analysis  of  this  data  ultimately 
requires  models  of  the  entire  set  of  cellular  pathways  on  a  quantitative  basis  which 
will  ultimately  reveal  important  feedback  mechanisms  that  are  undetectable  without 
this  comprehensive,  quantitative  approach.”  [25] 

Cell  models  must  be  incorporated  in  understanding  how  the  multiple  metabolic 
interactions  inside  a  cell  lead  to  changes  at  the  gene  and  protein  level.  Various 
research  papers  have  been  published  to  model  the  observed  functions  of  cells.  Over 
the  last  half  of  the  twentieth  century,  mathematics  has  also  played  a  large  role  in 
cellular  modeling.  “There  has  grown  a  sizable  literature  on  the  theoretical  analysis 
of  biological  regulation,  mostly  by  means  of  mathematical  models.”  [15]  As  more 
information  is  gained  about  the  internal  processes  of  cells,  mathematical  models  of 
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intracellular  reactions  have  been  incorporated  to  help  keep  track  of  these  cellular 
processes. 

Applications  of  mathematical  models  to  cell  processes  are  needed  for  many 
reasons.  They  are  needed  to  identify  the  design  principles  of  biochemical-based 
logic  within  a  cell.  This  is  key  in  understanding  how  normal  and  mutant  cells 
behave  in  response  to  external  and  internal  stimuli.  More  importantly,  they  can  be 
incorporated  into  simulations  that  will  give  an  understanding  of  how  cells  behave 
in  a  hypothesized  environment.  These  simulations,  combined  with  mathematical 
models,  will  predict  effects  of  different  nutrient  concentrations  or  cell  miitations  on 
regulatory  outcomes  of  the  cell.  This  will  lead  to  an  improved  understanding  of  cell 
toxicity  by  keeping  track  of  each  reaction  and  determining  how  the  cell  responds. 

The  cellular  processes  within  human  cells  have  great  sophistication  and  are 
extremely  complicated.  For  this  reason,  the  modeling  process  of  the  human  cell 
must  be  scaled  down  into  simplified  models  that  are  more  suitable  for  an  initial 
investigation. 

The  cell  is  comprised  of  many  complex  systems  that  control  cellular  functions, 
such  as,  metabolic  pathways,  transcription,  and  translation.  This  thesis  reviews 
work  previously  done  on  modeling  transcription,  translation,  and  metabolism  as  well 
as  theoretical  models  already  in  existence  such  as  the  Brusselator  and  Schnacken- 
berg’s  model.  Mathematically  modeling  small  systems  and  applying  them  to  larger 
cellular  processes  will  aid  in  understanding  how  cells  are  affected  by  toxic  chemicals. 

1.2  Problem 

In  order  for  the  Air  Force  to  develop  mathematical  models  to  analyze  cellular 
data,  it  is  important  to  understand  how  small  models  of  cellular  activity  were  created, 
as  well  as,  what  behavior  they  displayed  under  various  conditions.  By  constructing 
models  of  cellular  activity,  and  evaluating  their  behavior  in  response  to  a  sensitivity 
analysis,  larger  models  may  be  better  understood.  Understanding  the  underlying 
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behavior  of  small  models  and  their  relation  to  large  systems  will  lead  to  a  better 
understanding  of  how  the  Air  Force  should  construct  intracelhilar  models  to  assist 
in  future  toxicology  studies. 

1.3  Scope 

This  research  began  by  evaluating  cellular  models,  such  as  the  Hasty  model 
[14],  and  hjrpothetical  chemical  models,  such  as  the  Bnisselator  [27]  and  Schnacken- 
berg  [30]  models.  The  mathematical  models  presented  in  this  thesis  were  derived 
deterministically  from  reaction  equations  that  were  based  on  mass  action  kinetics. 
Each  reaction  equation  in  the  system  represents  complex  functions  of  molecular  in¬ 
teraction.  Reactants  within  each  equation  were  assumed  to  be  either  individiial 
molecules  or  groups  of  molecules.  The  rate  constants  within  the  reaction  equations 
reflect  an  average  rate  of  reaction.  The  models  are  considered  to  be  components  of 
larger  cellular  systems. 

1-4  Summary  of  Thesis 

This  thesis  is  organized  as  follows: 

Chapter  2  gives  a  brief  overview  of  cellular  components  and  functions.  It  then 
discusses  previous  and  current  work  on  the  development  of  mathematical  models 
of  cellular  activity.  It  introduces  computer  algorithms  used  to  describe  these  and 
larger  models.  It  also  introduces  hypothetical  chemical  oscillators  to  be  reviewed  as 
a  part  of  this  thesis  in  Chapters  3  and  4. 

Chapter  3  presents  techniques  for  describing  the  behavior  of  linear  and  nonlin¬ 
ear  systems.  It  then  derives  and  evaluates  the  behavior  of  current  models.  Finally, 
two  new  models  are  presented  as  examples  of  how  current  models  can  be  extended. 

Chapter  4  analyzes  the  behavioral  changes  of  current  and  new  models.  It 
then  addresses  how  some  types  of  behavior  may  cause  interesting  outcomes  within 
a  larger  cellular  model.  Finally,  a  small  model  is  incorporated  into  a  larger  model 
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in  an  effort  to  produce  a  greater  understanding  of  how  large  systems  are  affected  by 
small  component  models. 

Chapter  5  summarizes  the  work  completed.  It  also  presents  conclusions 
reached,  and  gives  recommendations  for  future  work. 


1-4 


II.  Background 


2. 1  OvervieAV 

Metabolic  pathways  can  be  described  both  qiialitatively  and  quantitatively. 
Experimental  observation  can  characterize  metabolic  pathways,  biit  a  mathematical 
approach  can  also  uncover  unknown  facts  about  the  process.  For  example,  mathe¬ 
matical  modeling  may  be  used  to  evaluate  concentration  levels  of  reactants  within  a 
chemical  or  biological  system  where  experiments  become  too  expensive,  dangerous, 
or  time  consuming  to  carry  out.  Together  these  two  types  of  descriptions  may  un¬ 
cover  unaccounted-for  or  misunderstood  processes  leading  to  a  better  understanding 
of  the  cell’s  processes.  Thus,  learning  about  cellular  processes  from  both  experi¬ 
mental  and  mathematical  modeling  is  a  very  beneficial  way  to  gain  understanding 
about  those  cellular  processes. 

Metabolic  regulation  is  an  important  cellular  process  typically  modeled  with 
nonlinear  ordinary  differential  equations  (ODE’s)  based  on  chemical  kinetics.  The 
need  for  metabolic  regulation  arises  as  a  result  of  different  stimuli  including  stress, 
changes  in  the  environment,  and  various  nutrients  ingested  into  the  cell.  Metabolic 
regulation  can  be  referred  to  as  signal  transmission  which  occur  through  chemical 
interactions,  enzymatic  reactions,  protein  degradation,  and  production  of  intracellu¬ 
lar  messengers.  It  addresses  questions  such  as  ‘How  do  pathways  change  to  increase 
production  of  products?’  and  ‘How  does  the  cell  respond  to  changes  in  nutrients?’ [20] 
These  pathways  are  linked  through  feedback  loops,  which  use  various  regulators  act¬ 
ing  as  messengers  that  are  produced  by  one  pathway  and  used  as  inputs  for  other 
pathways.  Ideally,  these  pathways  should  link  together  to  form  large  systems  of 
ODE’s  that  describe  the  cellular  process  as  a  whole.  “The  intricacy  and  variety  of 
biological  signaling  networks  often  defy  analysis  based  on  intuition.  Given  these  lui- 
certainties,  models  such  as  these  should  not  be  considered  as  definitive  descriptions 
of  networks  within  the  cell,  but  rather  as  one  approach  that  allows  us  to  understand 
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the  capabilities  of  complex  systems  and  devise  experiments  to  test  these  capabilities 

[4].” 


2.2  Biology 

Understanding  the  correct  approach  to  modeling  celhilar  behavior  results  from 
understanding  the  cellular  mechanism  being  modeled.  In  order  to  model  transcrip¬ 
tion  and  translation  of  a  prokaryotic  cell,  a  basic  understanding  of  the  biology  back¬ 
ground  behind  the  cell  must  be  outlined.  This  section  contains  qiiotes  that  address 
the  basic  definitions  for  eukaryotic  and  prokaryotic  organisms,  proteins  and  enzymes, 
DNA,  RNA  polymerase,  mRNA,  tRNA,  rRNA,  transcription,  translation,  operons, 
activators,  repressors,  and  metabolic  pathways  along  with  some  explanation  of  their 
functions.  Also  included  are  references  where  more  detail  can  be  found.  Once  these 
basic  definitions  are  addressed,  mathematical  models  of  some  of  these  biological 
systems  will  be  considered. 

2.2.1  Bacteriophage  A.  Bacteriophage  A  is  a  virus  infecting  E.  coli.  The 
bacteriophage  A  life  cycle  is  as  follows.  “The  infection  of  the  bacterial  host  E.  coli 
begins  when  the  virus  specifically  adsorbs  to  the  cell  and  injects  its  DNA.  The  lin¬ 
ear  DNA  then  circularizes  and  commences  directing  the  infection  process.  In  the 
lysogenic  mode,  the  phage  DNA  is  stably  integrated  at  a  specific  site  in  the  host 
chromosome  and  so  that  it  is  passively  replicated  with  the  bacterial  cell.  Alter¬ 
natively,  the  phage  may  take  up  the  lytic  mode  in  which  the  DNA  directs  its  own 
replication,  as  well  as  the  synthesis  of  viral  proteins  so  as  to  result  in  the  lysis  of  the 
host  cell  with  the  release  of  ~  100  progeny  phages.  DNA  damage,  as  is  caused, 
for  example,  by  UV  radiation,  induces  the  excision  of  the  prophage  DNA  from  the 
lysogenic  bacterial  chromosome  and  causes  the  phage  to  take  up  the  lytic  mode.” 
[34:p  1090]  More  information  about  bacteriophage  A  can  be  found  in  [34:p  1090] 
and  [19:pp  83,  92,  335-336]. 
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2.2.2  Prokaryotic  Organisms.  Prokaryotes  are  “organisms,  usually  bacte¬ 
ria,  that  have  neither  a  membrane-bound  nucleus  enclosing  their  chromosomes  nor 
functional  organelles  such  as  mitochondria  and  chloroplasts.”  [10:p  651]  More  in¬ 
formation  about  prokaryotic  organisms  can  be  foimd  in  [6:p  502],  [17:pp  25,  27,  52], 
[l:pp  20-21],  [10:ch  1,2],  and  [34], 

2.2.3  Eukaryotic  Organisms.  “Eukaryotes  are  defined  by  the  division  of 
each  cell  into  a  nucleus  that  contains  the  genetic  material,  surrounded  by  a  cyto¬ 
plasm,  which  in  turn  is  bounded  by  the  plasma  membrane  that  marks  the  periphery 
of  the  cell.  The  cytoplasm  contains  other  discrete  compartments,  also  bounded  by 
membranes.”  [17:p  27]  More  information  about  eukaryotic  organisms  can  be  found 
in  [6:p  6],  [l:p  237],  [10:ch  1,  2],  and  [34]. 

2.2.4-  Amino  Acids.  “Amino  acids  are  the  building  blocks  of  proteins.” 
[10;p  624]  “In  present-day  living  cells,  large  polypeptides-known  as  proteins-and 
polynucleotides-in  the  form  of  both  ribonucleic  acids  (RNA)  and  deoxyribonucleic 
acids  (DNA)-are  commonly  viewed  as  the  most  important  constituents.  A  restricted 
set  of  20  amino  acids  constitute  the  universal  building  blocks  of  the  proteins,  while 
RNA  and  DNA  molecules  are  constructed  from  just  four  types  of  nucleotides  each.” 
[l:p  4-5]  More  information  about  amino  acids  can  be  found  in  [6:pp  68-70],  [17:pp 
8-11],  [l:p  46],  [10:ch  1,  2],  and  [34]. 

2.2.5  Polypeptide.  A  polypeptide  is  “a  linear  series  of  amino  acids  linked 
together  with  peptide  bonds,  it  is  also  called  a  protein  or  protein  chain.”  [10:p  648] 
More  information  about  polypeptides  can  be  found  in  [6:pp  68-76],  [17:pp  3-25],  [l:pp 
107-108,  111-135],  [10:ch  1,2],  and  [34]. 

2.2.6  Metabolic  Pathways.  “As  a  whole,  metabolism  is  concerned  with 
managing  the  material  and  energy  resources  of  the  cell.  Some  metabolic  pathways 
release  energy  by  breaking  down  complex  molecules  to  simpler  compounds.  These 
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degenerative  processes  are  called  catabolic  pathways.  A  major  thoroughfare  of 
catabolism  is  cellular  respiration,  in  which  the  sugar  ghicose  and  other  organic  fuels 
are  broken  down  to  carbon  dioxide  and  water.  Energy  that  was  stored  in  the 
organic  molecules  becomes  available  to  do  the  work  of  the  cell.  Anabolic  pathways, 
in  contrast,  consume  energy  to  build  complicated  molecules  from  simpler  ones.  An 
example  of  anabolism  is  the  synthesis  of  a  protein  from  amino  acids.  Catabolic  and 
anabolic  pathways  are  the  downhill  and  uphill  avemies  of  the  metabolic  map.”  [6:pp 
83-84]  More  information  about  metabolic  pathways  can  be  found  in  [l:p  87],  [10:ch 
1,  2],  and  [34]. 

2.2. 7  DNA.  “Biological  instructions  are  encoded  in  the  molecules  known  as 
DN  A  (deoxyribonucleic  acid).  DNA  is  the  substance  of  genes,  the  units  of  inheritance 
that  transmit  information  from  parents  to  offspring.”  [6:pp  6,  77]  More  information 
about  DNA  can  be  found  in  [6:p  77],  [17:p  81],  [l:pp  4-5],  [10:ch  1,  2],  and  [34]. 

2.2.8  RNA  Polymerase.  “RNA  polymerase,  the  enzyme  responsible  for  the 
DNA-directed  synthesis  of  RNA,  was  discovered  independently  in  1960  by  Samual 
Weiss  and  Jerard  Hurwitz...All  cells  contain  RNA  polymerase.  In  bacteria,  one 
species  of  this  enzyme  synthesizes  all  of  the  cell’s  RNA  except  the  short  RNA  primers 
employed  in  DNA  replication.... RNA  polymerase,  which  has  a  characteristic  large 
size,  binds  to  DNA  as  a  protomer.  This  large  size  is  presumably  a  consequence  of 
the  haloenzyme’s  several  complex  functions  including  template  binding,  RNA  chain 
initiation,  chain  elongation,  and  chain  termination... Once  an  RNA  polymerase  mole¬ 
cule  has  initiated  transcription  and  moved  away  from  the  promoter,  another  RNA 
polymerase  can  follow  suit.  The  synthesis  of  RNAs  that  are  needed  in  large  quanti¬ 
ties,  ribosomal  RNAs,  for  example,  is  initiated  as  often  as  is  sterically  possible,  about 
once  per  second.”  [34:pp  919-920,  925]  More  information  about  RNA  polymerase 
can  be  found  in  [6:p  300],  [17:pp  811,  824],  [l:p  252],  [10:ch  1,  2],  and  [34]. 
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2.2.9  Transcription.  “Transcription  generates  a  single-stranded  RNA  iden¬ 
tical  in  sequence  with  one  of  the  strands  of  the  duplex  DNA.  Several  different  types 
of  RNA  are  generated  by  transcription;  the  three  principal  classes  involved  in  the 
synthesis  of  proteins  are:  messenger  RNA  (mRNA),  transfer  RNA  (tRNA),  and  ri- 
bosomal  RNA  (rRNA).”  [17:p  154]  “In  a  metabolically  active  cell,  about  3-5%  of 
the  cellular  RNA  is  mRNA,  about  90%  is  rRNA,  and  about  4%  is  tRNA.”  [10:p 
30]  Furthermore,  “the  in  vivo  rate  of  transcription  is  20  to  50  micleotides  per  sec¬ 
ond  at  37°  C  as  indicated  by  the  rate  at  E.  coli  incorporate  H-labeled  micleosides 
into  RNA.”  [34:p  925]  More  information  about  transcription  can  be  found  in  [6:pp 
296-297,  300-302,  315],  [l:p  366],  [10:ch  1,  2],  and  [34]. 

2.2.10  Operator.  The  operator  is  “the  region  of  DNA  that  is  upstream 
from  a  prokaryotic  gene  and  to  which  a  repressor  or  activator  binds.”  [10:p  647] 
More  information  about  operators  can  be  found  in  [6:p  338],  [17:p  166],  [l:p  417], 
[10:ch  1,  2],  and  [34]. 

2.2.11  Effector  Molecule.  An  effector  molecule  is  “a  low- molecular- weight 
compound  that  modifies  the  function  of  a  regulatory  protein.”  [10:p  634]  More 
information  about  effector  molecules  can  be  found  in  [10: ch  1,  2]. 

2.2.12  Activator.  An  activator  is:  “(1)  A  substance  or  physical  agent  that 
stimulates  transcription  of  a  specific  gene  or  operon.  (2)  A  protein  that  binds  to 
an  operator  and  enhances  the  rate  of  transcription.  Also  called  activator  protein.” 
[10:p  623]  More  information  about  activators  can  be  found  in  [6:pp  354-355],  [17:p 
359],  [10:ch  1,  2],  and  [34]. 

2.2.13  Activator  site.  An  activator  site  is  “a  DNA  sequence  to  which  an 
activator  protein  binds.  Also  called  activating  site.”  [10:p  623]  More  information 
about  activator  sites  can  be  found  in  [10:  ch  1,  2]. 
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2.2.14  Repressor.  A  repressor  is  “a  protein  that  binds  to  the  operator 
or  promoter  region  of  a  gene  and  prevents  transcription  by  blocking  the  binding  of 
RNA  polymerase.”  [10:p  653]  More  information  about  repressors  can  be  found  in 
[6:p  338],  [17:p  339],  [l:p  249],  [10:ch  1,  2],  and  [34], 

2.2.15  Operons.  “All  together,  the  operator,  the  promoter,  and  the  genes 
they  control-the  entire  stretch  of  DNA  required  for  enzyme  prodiiction  for  the  tryp¬ 
tophan  pathway-is  called  an  operon.”  [6:p  338]  More  information  about  operons 
can  be  found  in  [6:p  338],  [17:p  166],  [l:p  417],  [10:ch  1,  2],  and  [34], 

2.2.16  Messenger  RNA.  “Each  gene  along  the  length  of  the  DNA  molecule 
directs  the  synthesis  of  a  type  of  RNA  called  messenger  RNA  (mRNA).  The  mRNA 
molecule  then  interacts  with  the  cell’s  protein-synthesizing  machinery  to  direct  the 
production  of  a  polypeptide.  We  can  summarize  the  flow  of  genetic  information  as 
DNA^RNA^protein.”  [6:p  77]  More  information  about  messenger  RNA  can  be 
found  in  [6:p  296],  [17:p  717],  [l:p  223],  [10:ch  1,  2],  and  [34]. 

2.2.17  Ribosomal  RNA.  “There  are  two  major  types  of  rRNA.  The  larger 
of  these  rRNA’s  combines  with  a  set  of  proteins  to  form  a  ribonucleoprotein  complex 
called  the  large  ribosomal  subunit.  The  smaller  rRNA  combines  with  another  set  of 
proteins  to  form  a  smaller  ribosomal  subunit.  During  protein  synthesis,  one  large 
ribosomal  subunit  and  one  small  ribosomal  subunit  combine  to  form  a  ribosome.” 
[10:p  30]  Furthermore,  “even  E.  coli  needs  seven  copies  of  its  rRNA  genes  to  keep 
up  with  the  cell’s  need  for  ribosomes.”  [l:p  378]  More  information  about  ribosomal 
RNA  can  be  found  in  [6:p  306],  [17;p  179],  [lOxh  1,  2],  and  [34]. 

2.2.18  Ribosome.  “The  actual  sites  of  protein  synthesis  are  cellular  struc¬ 
tures  called  ribosomes.”  [6:p  77]  More  information  about  Ribosomes  can  be  found 
in  [6:pp  111,  306,  506,  523],  [17:pp  6,  33,  159-162],  [l:p  223],  [10:ch  2],  and  [34]. 
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2.2.19  Transfer  RNA.  “In  the  process  of  translation  a  cell  interprets 
a  genetic  message  and  builds  a  protein  accordingly.  The  message  is  a  series  of 
codons  along  an  mRNA  molecule,  and  the  interpreter  is  transfer  RNA  (tRNA).  The 
function  of  tRNA  is  to  transfer  amino  acids  from  the  cytoplasm’s  amino  acid  pool 
to  a  ribosome.”  [6:p  304]  More  information  aboiit  transfer  RNA  can  be  found  in 
[17:pp  153,  213],  [l:pp  227-230],  [lOxh  1,  2],  and  [34]. 

2.2.20  Translation.  “Translation  converts  the  nucleotide  sequence  of  RNA 
into  the  sequence  of  amino  acids  comprising  a  protein.  An  mRNA  is  translated  into 
a  protein  sequence:  tRNA  and  rRNA  provide  other  components  of  the  apparatus 
for  protein  synthesis.  The  entire  length  of  an  mRNA  is  not  translated,  but  each 
mRNA  contains  at  least  one  coding  region  that  is  related  to  a  protein  seqiience  by 
the  genetic  code:  each  nucleotide  triplet  (codon)  of  the  coding  region  represents  one 
amino  acid.”  [17:p  154]  More  information  about  translation  can  be  found  in  [6:pp 
296-297,  304],  [l:p  25],  [10:ch  1,  2],  and  [34]. 

2.3  E-Cell 

The  E-Cell  project  began  in  1996  at  Keio  University  [33].  It’s  goal  is  to  model 
various  biochemical  processes  with  the  intended  goal  of  eventually  modeling  the 
entire  cell.  E-Cell’s  main  attraction  is  its  ability  to  integrate  metabolic  pathways 
with  high-order  cellular  processes.  These  processes  include  protein  synthesis  through 
transcription  and  translation  and  membrane  transport. 

E-Cell  is  run  under  a  UNIX,  Linux,  MSDOS,  or  MS  Windows  operating  system. 
It  requires  the  user  to  define  the  cell’s  molecules,  locations,  and  concentrations.  It 
then  requires  the  user  to  input  the  set  of  reaction  rules  that  govern  all  these  processes. 
E-Cell  then  computes  the  concentration  and  location  of  each  molecule  at  each  time 
increment.  Eurt  her  more,  in  UNIX  or  Linux,  a  graphical  interface  allows  the  user  to 
monitor  the  cell  as  each  reaction  occurs. 
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Using  E-Cell,  Tomita[33]  et  al.  were  able  to  construct  an  electronic  cell  with  127 
genes  that  was  capable  of  “self-support.”  The  cell  model  included  “genes  for  tran¬ 
scription,  translation,  glycolysis  pathways  for  energy  prodiiction,  membrane  trans¬ 
port,  and  the  phospholipid  biosynthesis  pathway  for  membrane  structure.”  [33] 

2.4  McAdams,  Arkin,  Shapiro,  and  Bhalla 

Biochemical  and  genetic  approaches  have  identified  a  mimber  of  the  molecular 
mechanisms  of  bacteria.  A  problem  attracting  considerable  research  effort  is  the 
detailed  understanding  of  the  dynamic  behavior  of  large  systems  of  genes  and  related 
proteins  and  how  they  interact  to  control  cellular  fimctions  over  the  cell’s  life  cycle. 
Harley  McAdams  and  Lucy  Shapiro  [22]  investigated  this  problem  through  the  use 
of  genetic  networks.  These  genetic  networks  attempt  to  incorporate  some  of  the 
techniques  used  in  electrical  circuits  to  gain  a  better  understanding  of  intracellular 
mechanisms.  Genetic  networks  consist  of  hundreds  of  genes  and  consequently  are 
complex.  “Extension  of  metabolic  modeling  methods  to  include  more  realistic  ge¬ 
netic  regulatory  mechanisms  is  a  current  challenge  to  the  field.”  [20]  For  instance,  the 
coupled  reactions  controlling  cell  division  in  prokaryotic  cells  have  not  been  identi¬ 
fied.  Due  to  the  commonalities  in  the  function  of  these  biochemically  based  genetic 
circuits  and  electrical  circuits  [20] ,  a  new  modeling  technique  has  emerged  as  a  way 
of  integrating  conventional  biochemical  kinetic  modeling  within  the  framework  of  a 
circuit  simulation. 

In  genetic  networks,  the  protein  production  encoded  by  one  gene  often  regulates 
expression  of  other  genes.  Activation  depends  on  the  level  of  expression  of  one 
gene  required  to  control  the  expression  of  other  genes.  The  concept  of  genetic 
circuits  was  motivated  by  electrical  circuits,  but  they  have  important  differences. 
McAdams  and  Shapiro  state  that,  “Common  transistor  circuits  can  operate  at  more 
than  10®  cycles  per  second.  In  contrast,  the  protein  signal-controlled  switching  rate 
in  genetic  circuits  is  around  10^^  per  second.”  [22]  Due  to  this,  McAdams  and 
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Shapiro  incorporated  time  delays  into  their  system.  They  preformed  simulations 
of  gene  expression  which  resulted  in  the  proteins  in  the  simiilation  being  produced 
in  short  bursts  and  in  variable  numbers  occurring  at  random  times.  As  a  result, 
they  claim  there  can  be  large  gaps  in  time  between  siiccessive  events  in  regulatory 
changes  across  a  cell.  Additionally,  randomized  patterns  of  expression  of  competitive 
effector  molecules  can  produce  outcomes  in  switching  mechanisms  that  select  between 
alternative  regulatory  paths.  As  a  result  of  following  different  paths,  different  types 
of  cells  may  occur.  McAdams  and  Shapiro  observed  this  in  their  model  of  the 
Bacteriophage  A  lysis-lysogeny  decision  circuit. 

Molecules  that  control  gene  expression,  such  as  activators  and  repressors,  may 
act  at  extremely  low  intracellular  concentrations.  For  this  reason,  large  fluctuations 
in  reaction  rates  may  occur,  causing  variations  of  concentrations  of  each  molecular 
species.  In  spite  of  this,  many  regulatory  pathways  in  cells  have  highly  predictable 
outcomes.  To  achieve  regulatory  reliability,  cells  use  redundant  genes  and  networks, 
as  well  as,  feedback  loops  to  buffer  pathways  against  mutational  or  environmental 
perturbations.  [21] 

Bhalla  and  Iyengar  [4]  have  also  studied  metabolic  pathways  using  circuit  com¬ 
ponents.  They  have  shown  that  many  distinct  signaling  pathways  allow  the  cell  to 
receive  information,  process  it,  and  respond.  Their  results  include  networks  that  ex¬ 
hibit  properties  such  as  integration  of  signals  across  multiple  time  scales,  generation 
of  distinct  outputs  depending  on  input  strength  and  duration,  and  self-sustaining 
feedback  loops.  They  state  that,  “feedback  loops  can  result  in  bistable  behav¬ 
ior  with  discrete  steady-state  activities,  well-defined  input  thresholds  for  transition 
between  states  and  prolonged  signal  output,  and  signal  modulation  in  response  to 
transient  stimuli.”  [4] 
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2.5  Smolen  et  al.  and  Hasty  et  al. 

2.5.1  Smolen  et  al.  Smolen  et  al.  [32]  investigate  the  possibility  that 
multi-stability  and  oscillatory  behavior  occur  in  genetic  regulatory  systems.  They 
created  kinetic  models  that  incorporate  known  features  of  genetic  regulatory  sys¬ 
tems.  “These  include  autoregulation  and  stimulus-dependent  phosphorylation  of 
transcription  factors  (TFs),  dimerization  of  TFs,  crosstalk,  and  feedback.  The  sim¬ 
plest  model  manifested  multiple  stable  steady  states,  and  brief  perturbations  could 
switch  the  model  between  these  states.  ...  In  slightly  more  complex  models,  oscilla¬ 
tory  regimes  were  identified.”  [32] 

Smolen  et  al.  looked  at  different  complexities  of  models  to  determine  where 
strange  behavior  in  genetic  systems  may  occur.  They  began  by  considering  “a  rel¬ 
atively  simple  model  of  transcription  factors  (TFs)  subject  to  positive  and  negative 
autoregulation  of  their  own  transcription.”  [32]  First,  they  considered  “signal- 
transduction  pathways  in  which  stimuli  lead  to  second  messenger  generation  and 
phosphorylation  of  TF’s,  which  in  turn  bind  to  DNA  sequences  known  as  responsive 
elements  and  thereby  regulate  the  transcription  of  specific  genes.  The  regulatory 
activity  of  TFs  is  often  modulated  by  phosphorylation  and  by  intermolecular  inter¬ 
actions.  For  example,  TFs  often  bind  to  DNA  as  homodimers  or  as  hetrodimers 
of  different  TF  family  members.. .As  a  result,  some  TFs,  such  as  Jun,  autoregulate 
their  own  transcription  [24]... Thus  responsive  elements  affecting  TF  gene  transcrip¬ 
tion  can  provide  crosstalk  and  positive  feedback.”  [32]  Smolen  et  al.  made  some 
simplifications  to  the  model  so  it  could  be  appreciated  intuitively  and  described 
mathematically.  “A  single  transcriptional  activator,  which  we  term  TF-A,  is  con¬ 
sidered  as  part  of  a  pathway  mediating  a  cellular  response  to  a  stimulus.  The  TF 
forms  a  homodimer  that  can  bind  to  responsive  elements  (TF-REs).  The  tf-a  gene 
incorporate  one  of  these  responsive  elements,  and  when  homodimers  bind  to  this 
element  TF-A  transcription  is  increased.  Binding  to  the  TF-REs  is  independent  of 
dimer  phosphorylation.  Only  phosphorylated  dimers,  however,  can  activate  tran- 


2-10 


scription.  The  fraction  of  dimers  phosphorylated  is  dependent  on  the  activity  of 
kinases  and  phosphatases  whose  activity  can  be  regulated  by  external  signals.  Thus 
this  model  incorporate  both  signal-activated  transcription  and  positive  feedback  on 
the  rate  of  TF  synthesis.”  [32] 

Smolen  et  al.  expanded  their  model  by  incorporating  both  positive  and  nega¬ 
tive  feedback.  “Responsive  elements  have  also  been  found  that  regulate  genes  for 
potent  transcriptional  inhibitors... Such  responsive  elements  provide  negative  feed¬ 
back  loops.  Given  the  above  interactions,  there  is  the  possibility  for  rich  dynamic 
activity. .However,  models  similar  to  the  one  above  have  only  one  type  of  feedback: 
either  positive  or  negative.  Without  both  types  of  feedback,  siich  a  model  cannot 
be  expected  to  support  oscillations.  For  example,  positive  feedback  can  act  to  drive 
[TF-A]  to  high  levels,  but  then  there  is  no  process  to  bring  [TF-A]  back  down.  To 
investigate  the  effects  of  negative  feedback,  we  introduced  a  protein,  TF-R,  that 
represses  transcription  by  binding  to  the  TF-REs.  Its  rate  of  synthesis  is  increased 
by  binding  of  the  TF-A  dimer  to  a  TF-RE.  TF-R  competitively  inhibits  binding 
of  TF-A  dimers  to  TF-REs;  thus  it  inhibits  the  transcription  of  the  genes  tf-a  and 
tf-r.”  [32] 

Using  a  bifurcation  graph  that  plotted  steady  state  concentrations  of  the  tran¬ 
scription  factor  as  a  function  of  a  parameter,  Smolen  et  al.  determined  that  de¬ 
pending  on  different  parameters,  the  steady  state  concentrations  of  this  model  may 
exhibit  strange  behavior  in  the  form  of  two  stable  steady  states  and  one  unstable 
steady  state.  This  behavior  is  known  as  bistability.  Smolen  et  al.  stated  that 
bistable  “transitions  could  correspond  physiologically  to  brief  stimuli,  such  as  ex¬ 
posure  to  hormone,  leading  to  long-lasting  increases  or  decreases  in  the  levels  of 
particular  proteins.”  [32] 

2.5.2  Hasty  et  al.  Hasty  et  al.[14]  based  their  work,  in  part,  on  that  of 
Smolen  et  al.  [32]  In  their  analysis,  the  original  model  created  by  Smolen  et  al. 
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that  included  only  the  activator  was  incorporated  into  the  following  system  which 
resembles  the  stoichiometric  equations  of  chemistry. 


2X  ^  X2  (2.1) 

fc_l 

D  +  X2^  DX2  (2.2) 

_ 2 

DX2  +  P  DX2  +  nX  (2.3) 

X  A  (2.4) 


In  the  presentation  by  Hasty  et  ah,  they  assume  the  following.  “Let  X  denote 
the  protein  and  D  a  specific  binding  site  on  a  DNA  promoter.  The  prodiiction  of  X 
takes  place  when  the  dimer  X2  forms  a  complex  with  the  DNA  at  site  D...The  first 
reaction  represents  the  dimerization  of  X,  and  the  second  the  formation  of  a  DNA- 
promoter  complex.  In  the  presence  of  RNA  polymerase  P  and  the  reactants  required 
for  transcription,  this  complex  leads  to  the  effective  transcription  of  the  gene  coding 
for  X.  The  produced  mRNA  is  then  translated  into  protein  X,  where  n  is  the  number 
of  proteins  per  transcript.  Both  transcription  and  translation  are  represented  by  the 
third  equation.  Lastly  X  degrades  via  reaction  four... Further,  the  first  two  reactions 
are  orders  of  magnitude  faster  than  transcription  and  degradation.”  [14] 

Analysis  led  Hasty  et  al.  to  conclude  that  “bistability  arises  as  a  consequence 
of  both  the  dimerization  of  the  protein,  and  the  competition  between  the  production 
via  transcription  and  the  degradation  in  the  last  reaction.”  [14] 

In  the  second  paper  by  Hasty  et  al.[13]  they  model  the  effects  the  bacteriophage 
A  virus  by  analyzing  the  evolution  of  the  A  repressor  protein.  “Bacteriophage  A  has 
two  alternate  life  cycles,  lytic  and  lysogenic  growth.”  [19:83,  92,  335-336]  “In  the 
context  of  the  lysis-lysogeny  pathway  in  the  A  virus,  the  autoregulation  of  A  repressor 
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expression  is  well  characterized  [16],  [28]...  We  present  two  models  describing  the 
regulation  of  such  a  network.”  [13] 

Hasty  et  al.  constructed  the  model  to  represent  the  different  effects  bistability 
has  on  mutant  and  normal  operator  regions  of  the  bacteriophage  A  virus.  “Although 
the  full  promoter  region  in  A  phage  contains  the  three  operator  sites  known  as  ORl, 
OR2,  and  OR3,  we  first  consider  a  mutant  system  whereby  the  operator  site  ORl 
is  absent  from  the  region.  The  basic  dynamical  properties  of  this  network,  along 
with  a  categorization  of  the  biochemical  reactions,  are  as  follows  [16],  [28].  The 
gene  cl  expresses  repressor  (Cl),  which  in  turn  dimerizes  and  binds  to  the  DNA  as  a 
transcription  factor.  In  the  mutant  system,  this  binding  can  take  place  at  one  of  the 
two  binding  sites,  OR2  or  OR3.  (Here  we  ignore  nonspecific  binding.)  Binding  at 
OR2  enhances  transcription,  which  takes  place  downstream  of  OR3,  whereas  binding 
at  OR3  represses  transcription,  effectively  turning  off  production.”  [13] 

The  reactions  describing  the  mutant  operator  region  model  are  below. 


kl 

2X^X2 

(2.5) 

k2 

D  +  X2-  DX2 

(2.6) 

D  +  A2  S  DX; 

(2.7) 

DX2  +  A2  S  DX2X2 

(2.8) 

DX2  +  P^  DX2  +  P  +  nX 

(2.9) 

A  “  A 

(2.10) 

Equation  (2.5)  of  the  new  model  involves  the  dimerization  of  the  A 

protein  (X)  to  create  a  dimer  protein.  The  dimer  protein  {X2)  can  bind  to 

(D)  at  one  of  three  operator  sites.  The  second  operator  site,  labeled  as 

repressor 

the  DNA 

{DX2)  in 

the  model  above,  enhances  transcription,  which  is  consistent  with  lysogenic  growth 
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of  Bacteriophage  A.  The  third  operator  site,  labeled  as  [DX^)  in  the  model  above, 
which  exists  down  stream  from  the  second  operator  site  represses  transcription,  which 
is  consistent  with  lytic  growth  of  Bacteriophage  A.  Also,  a  dimer  protein  can  bond 
to  the  third  site,  and  then  the  second  {DX2X2).  In  this  case,  the  second  site 
acts  as  a  repressor  and  blocks  the  activator  site  from  moving  on.  Finally,  RNA 
polymerase  joins  the  activated  site  to  complete  transcription,  after  which,  translation 
is  completed,  and  as  a  last  process,  the  repressor  protein  degrades.  These  are  much 
slower  reactions  that  the  first  four  and  are  also  irreversible. 

As  in  the  previous  model.  Hasty  et  al.  states  that  the  bistability  in  this  system 
arises  from  the  competition  between  the  production  of  the  repressor  protein  seen  in 
equation  (2.9),  its  dimerization  in  equation  (2.5),  and  its  degradation  from  eqiiation 

(2.10). 

2.6  Chen  et  al. 

Chen,  He,  and  Church[7]  formulated  a  different  approach  to  studying  large 
systems  of  genes.  They  explain  how  to  construct  their  model  from  initial  data  given 
by  samples  of  mRNAs  and  proteins.  Their  results  suggest  that  given  accurate  data, 
their  models  can  determine  most  gene  regulation  accurately  at  the  genome  level. 

They  begin  by  considering  the  simplifications  that  must  be  made  about  the 
actual  processes  of  transcription  and  translation  in  order  to  fit  them  into  the  model. 
They  ignore  direct  feedback  from  mRNAs  to  genes.  Instead  they  assumed  all  feed¬ 
back  to  the  gene  is  a  direct  result  of  proteins.  Furthermore,  they  assumed  “the 
translation  mechanism  is  relatively  stable  (at  least  for  a  short  time),  so  the  feedback 
from  proteins  to  mRNAs  has  no  effect.  Each  mRNA  and  protein  molecule  degrades 
randomly,  and  its  components  are  recycled  in  the  cell.”  [7]  Finally,  they  assume 
there  is  no  feedback  from  the  metabolic  pathways  to  transcription.  They  point  out 
that  feedback  from  the  metabolic  pathway  to  transcription  may  actually  have  some 
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effect  on  the  system.  However,  it  is  a  simplification  they  have  chosen  to  more  easily 
solve  their  model. 

After  making  their  assumptions,  they  were  able  to  mathematically  describe 
their  model  through  the  nonlinear  dynamic  system  below. 


dr 

dt 

dp 

dt 


/(p)  -  Vr 
Lr-Up 


(2.11) 


The  variables  r  and  p  are  functions  of  time  t.  The  quantity  n  is  the  number 
of  genes  being  included  in  (2.11).  r  is  the  n-dimensional  vector  valued  function  for 
mRNA  concentrations,  p  is  the  n-dimensional  vector  valued  function  for  protein 
concentrations.  /  is  the  set  of  functions  of  proteins  p  that  describe  transcription. 
L  is  the  n  x  n  dimensional  matrix  describing  translation.  V  is  the  n  x  n  dimensional 
matrix  describing  degradation  of  mRNAs.  Finally,  U  is  the  n  x  n  dimensional  matrix 
describing  degradation  of  proteins. 

For  equations  (2.11),  the  following  observations  can  be  made.  mRNA  concen¬ 
trations  {r}  increase  as  transcription  {/(p)}  increases,  and  decreases  as  the  degra¬ 
dation  of  mRNA  {Vr}  increases.  Protein  concentrations  {p}  increase  as  translation 
of  mRNA  {Lr}  increases  and  decreases  as  degradation  of  proteins  {Up}  increases. 
Translation  rates  and  degradation  rates  are  assumed  constant  for  both  mRNA  and 
proteins  so  L,  V,  and  U  are  diagonal  matrices  with  no  zeros  on  the  diagonal.  Finally, 
they  assume  that  for  this  model,  there  is  no  time  delay.  They  leave  time  delay  for 
another  model. 
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They  solve  their  problem  by  linearizing  the  transcription  functions,  /(p).  The 
resulting  linear  system  has  a  unique  solution  which  can  be  written  in  terms  of  the 
eigenvalues  and  eigenvectors  of  the  system. 

Additional  insight  concerning  this  problem  may  be  found  by  considering  gene 
regulatory  networks.  Although  genes  act  as  regulators  for  other  genes,  this  regula¬ 
tion  often  involves  only  a  few  number  of  regulators.  Hence,  all  genes  do  not  act  as 
regulators  for  each  other.  Chen  claims  “that  the  number  of  regulators  for  a  gene 
is  small,  mostly  less  than  10  [7].”  Savageau  adds  validity  to  this  claim  by  stating 
“molecular  analysis  of  gene  regulation  in  bacteria  has  shown  that  most  gene  circuits 
are  governed  by  a  small  number  of  regulators  usually  one  to  three.  In  eukaryotes  the 
numbers  are  larger  in  some  cases,  but  seldom  more  than  a  dozen  regulators  influence 
a  given  gene  circuit.... The  most  biologically- suggestive  behaviors  were  found  when 
each  circuit  was  subject  to  two  or  three  regulatory  interactions  [29].”  This  fact 
allows  for  the  production  of  a  sparse  matrix. 

In  discussing  the  limitations  of  their  model  (2.11),  Chen  et  al.  make  note  of 
speciflc  potential  flaws  in  the  model.  They  point  out  that  although  the  model  does 
not  include  time  delay,  it  does  capture  many  features  of  gene  expression.  Further¬ 
more,  the  main  design  flaw  comes  from  the  ignorance  of  other  regulators  such  as 
those  constructed  via  the  metabolic  pathway.  Other  potential  hazards  of  this  model 
include  the  assumption  of  periodic  cell  cycles  and  small  numbers  of  regulators.  Fail¬ 
ure  to  meet  both  assumptions  prevents  the  model  from  being  reconstructed  via  the 
data. 

2.1  Heinrich  et  al. 

Heinrich  et  al.  [15]  described  metabolic  systems  through  mathematics  and  chem¬ 
ical  kinetics.  They  concluded  that  the  number  of  hypothetical  biological  concepts 
required  to  explain  data  results  may  be  scaled  down  by  the  use  of  models.  They 
explained  how  to  set  up  and  simplify  a  model  by  only  including  key  elements.  Then, 


2-16 


how  to  construct  the  model,  using  equations  based  on  thermodynamics  and  kinetic 
equations.  They  analyzed  the  dynamics  behind  these  models,  and  identified  princi¬ 
ples  of  stability.  Their  rational  for  this  analysis  was  to  show  how  biological  dynamics 
can  be  better  understood  when  incorporated  with  mathematical  models. 

Heinrich  et  al.  began  by  describing  the  basic  need  for  mathematical  modeling 
in  the  metabolic  pathway.  Modeling  can  be  used  to  discover  the  fate  of  some  of  the 
intermediary  proteins  and  enzymes  in  the  metabolic  pathway.  These  proteins  and 
enzymes  may  behave  differently  under  different  concentrations.  To  understand  these 
differences,  models  are  needed  to  analyze  the  basic  dynamics  and  stability  regions 
of  these  intermediary  protein  and  enzyme  concentrations. 

To  deal  with  the  complexity  of  biological  entities,  Heinrich  et  al.  describe  how 
to  simplify  models  based  on  the  models’  processes.  The  rate  at  which  regulatory 
properties  of  metabolism  are  expressed  depend,  in  part,  on  the  components  within 
metabolism,  concentration  of  the  different  elements,  and  the  kinetics  governing  the 
reactions.  It  is  important  to  understand  which  components  are  necessary  to  include 
in  the  model,  and  how  they  are  connected.  Heinrich  et  al.  explained  that  “only 
the  essential  reacting  variables  have  to  be  considered  since  faster  variables  are  of¬ 
ten  in  a  quasi-steady  state  and  slower  ones  are  constant.  ...  Whole  pathways  may 
be  substituted  by  single  reactions  in  this  way.  For  instance,  Selkov  [31]  considered 
both  hexokinase  and  phosphofructokinase  as  a  single  phosphorylating  reaction  and 
the  phosphoglycerate  kinase  and  pyruvate  kinase  reaction  as  one  reaction  regener¬ 
ating  ATP.”  [15]  Many  reactions  are  quite  complicated,  especially  when  enzymes 
are  involved.  Reactions  may  fluctuate  either  quickly  or  slowly  depending  on  the 
concentration  of  the  reactants.  For  this  reason,  simplified  rate  functions  are  often 
needed  to  describe  the  effects  of  the  different  concentrations  of  reactants.  They  are 
incorporated  into  systems  of  first  order  differential  equations  that  represent  the  flux 
in  concentrations  of  all  the  reactants  being  modeled.  This  flux  is  determined  by  the 
equations  below. 
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(2.12) 


^  CijVj  where  i  =  (1, n) 

1=1 

The  variables  in  equation  (2.12)  are  describes  as  follows:  Si  is  the  concentration 
of  the  reactant;  n  is  the  number  of  reactants;  r  is  the  mimber  of  reactions  in  the 
model;  Vj  is  the  combination  of  different  reactants  associated  with  S'*;  is  the 
number  of  molecules  of  reactant  Si  participating  in  reaction  j;  and  t  is  the  time. 
Now,  Cij  is  positive  if  reactant  Si  is  a  substrate  of  reaction  j,  where  as  Cij  is  negative 
if  reactant  Si  is  a  product  of  reaction  j.  Furthermore,  vj  can  be  determined  via  the 
following  equation. 

Vj  =  vJ  -  =  {kj  JJ  S-’^  -  JJ  S-’^)Rj{Si,pk)  (2.13) 

i  i 

products  substrates 

Equation  (2.13),  is  the  forward  reaction  rate  and  vJ  is  the  reverse  reaction 
rate.  Also,  kj  is  the  forward  reaction  rate  coefficient  and  k^j  is  the  reverse  reaction 
rate  coefficient.  Rj  is  the  function  containing  all  the  information  related  to  special 
catalytic  rate  functions  such  as  the  Michaelis- Menton  equation.  Lastly,  pk  is  the 
kinetic  parameters  associated  with  special  catalytic  rate  functions  such  as  iLm- values. 

Equation  (2.12)  can  be  represented  in  the  following  abbreviated  manner, 

=  fi{Si, ...,  Sn)  where  {i  =  1, ...,  n)  (2.14) 

where  fi  is  the  flux  of  reactant  Si  for  {i  =  1, ...,  n). 

An  analysis  of  this  mathematical  model  can  be  performed  to  determine  the 
dynamics  behind  the  model,  and  identify  principles  of  stability.  The  steady  state  of 
system  (2.14)  is  achieved  when  ^  =  0  =  fi{Sf, ...,  S')(),  where  S'^defines  the  steady 
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state  concentrations  of  the  individual  element  Si.  A  discussion  of  stability  will  be 
presented  in  Chapter  3  when  specific  applications  are  considered. 

The  domain  surrounding  a  steady  state  point  may  be  stable,  or  unstable.  If 
a  concentration  is  pulled  away  from  a  stable  steady  state,  it  will  eventually  move 
back  to  that  steady  state.  If  the  concentration  is  pulled  away  from  an  unstable 
state,  it  continues  to  move  away.  Although  biological  organisms  generally  involve 
stable  domains,  in  nonlinear  systems,  locally  stable  domains  may  be  very  small. 
Thus,  a  small  perturbation  from  the  stable  state  may  allow  it  to  return,  biit  a  larger 
perturbation  may  move  it  out  of  the  local  domain  and  in  the  direction  of  another 
steady  state,  a  limit  cycle,  or  cause  it  to  diverge. 

Small  systems  that  can  be  evaluated  have  the  capability  to  be  analyzed  via 
bifurcation  diagrams  in  the  form  of  plots  that  represent  steady  state  concentration 
as  a  function  of  a  parameter.  If  the  parameter  is  changed,  then  the  steady  state 
concentrations  may  increase  or  decrease.  Given  enough  change  in  the  parameter,  the 
steady  state  concentration  may  fork  into  multiple  steady  states  such  as  a  model  by 
Baras  et  al.  [3]  shows,  or  become  oscillatory  such  as  Gillespie  [9]  point  out  though  use 
of  the  Brusselator  equation.  This  diagram  is  critical  in  analyzing  dynamic  systems 
because  it  can  display  changes  of  the  stability  of  the  system  as  parameters  change. 

To  prove  their  case,  Heinrich  et  al.  [15]  evaluated  a  small  system  without 
linearizing  it.  They  obtained  numerical  solutions  of  the  system.  The  system 
exhibited  interesting  behavior  such  as  limit  cycles  and  multiple  stable  states.  This 
system  was  an  approximate  representation  of  a  form  of  metabolic  control  called 
feedback  loops.  One  variation  of  a  feedback  loop  is  where  a  product  of  a  pathway 
is  used  as  a  rate  controlling  mechanism  to  control  its  own  pathway.  This  type 
of  feedback  can  cause  the  pathway  to  have  multiple  steady  states.  Thus,  changes 
in  concentrations  may  lead  to  different  steady  states  and  eventually  to  different 
biological  behavior. 
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2.8  The  Brusselator  a,nd  Schnakenberg  equations 

The  Brusselator  model  was  created  as  “a  scheme  which  may  lead  to  a  chemical 
symmetry-breaking  instability.”  [27]  It  involves  four  reaction  eqiiations  with  two 
intermediate  reacting  components.  One  of  the  components  acts  as  an  autocatalyst. 
That  is,  it  catalyzes  its  own  formation.  In  doing  so,  it  involves  a  tri-molecular 
reaction  step.  Initially  introduced  as  a  chemical  oscillator,  Prigogine  and  Lefever 
considered  it  “physically  unrealistic  because  of  the  tri-molecular  step.”  [27]  How¬ 
ever,  it  provides  a  basic  understanding  of  how  oscillation  occurs,  and  may  lead  to 
a  greater  understanding  of  oscillations  in  biological  systems  when  incorporated  into 
more  general  systems  of  reaction  equations. 

Schnakenberg  [30]  considered  an  extremely  simple  mathematical  model  display¬ 
ing  oscillatory  limit  cycle  behavior.  He  began  his  analysis  by  considering  a  model  of 
general  chemical  reactions.  Systems  of  chemical  reactions  consist  of  stoichiometric 
equations.  Stoichiometric  systems  of  equations  take  the  form 


Aj  -|-  cijiS*!  -|-  C2jS2  +  csjS's  -f  ...  ^  Bj  -f  C2jS2  +  +  ...  (2.15) 

where  A  and  B  are  fixed  concentrations  of  external  reactants,  and  all  other  variables 
are  described  in  equation  (2.12) .  To  obtain  the  desired  mix  of  simplicity  in  the  model 
and  interesting  behavior  in  the  solution,  Schnackenberg  chose  to  use  two  variables 
of  reactants  to  create  his  model  which  is  of  the  form  of  equation  (2.14)  with  n  =  2. 
Schnackenberg  linearized  his  nonlinear  system  and  solved  to  find  the  eigenvalues  of 
the  linearized  system.  This  analysis  will  be  presented  in  Chapter  3  when  several 
examples  are  considered. 

With  the  understanding  that  at  least  one  unstable  steady  state  point  is  needed 
to  get  oscillatory  behavior,  Schnackenberg  continued  the  process  of  creating  a  simple 
oscillator.  He  did  this  by  first  limiting  the  number  of  equations  he  uses.  Systems 
of  four  reactions,  such  as  the  Brussilator,  have  already  been  modeled,  so  to  model 
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a  system  with  more  reactions  than  three  wouldn’t  lead  to  a  simplest  oscillatory 
equation.  Systems  of  only  a  single  reaction  can  be  ruled  out  because  their  steady 
state  is  the  thermodynamic  equilibrium.  For  this  reason,  oscillation  cannot  occur. 
Systems  with  two  equations  and  two  variables  can  also  be  ruled  out  because  they 
are  purely  stable  or  unstable.  Unstable  steady  states  in  this  system  will  only  create 
harmonic  instability  that  degenerates  out  to  infinity.  Thiis,  a  three  reaction  system 
is  considered. 

All  three  equations  can  be  combined  to  produce  the  reqiiired  conditions  of  a 
unstable  steady  state.  Thus,  completing  conditions  to  prodiice  a  limit  cycle  with 
only  three  equations  and  two  variables.  In  producing  a  limit  cycle,  it  is  important 
to  be  aware  that  the  conditions  applied  to  create  a  limit  cycle  don’t  rule  out  the 
possibility  that  both  an  unstable  and  a  stable  steady  state  lie  within  or  outside  the 
limit  cycle.  In  some  cases,  all  the  conditions  for  limit  cycles  may  exist,  but  if  there  is 
a  stable  node  inside  this  region  as  well,  the  concentrations  may  approach  that  stable 
node  instead  of  the  limit  cycle. 

Limit  cycles  are  of  great  importance  in  biology  and  particularly  in  metabolic 
pathways.  They  were  discussed  back  in  1968  by  E.  Selkov  [31]  as  models  for  self¬ 
oscillations  in  glycolysis.  It  is  not  unreasonable  to  think  that  such  models  may  be 
extended  to  other  systems  in  biology  to  include  transcription  and  translation. 

2. 9  Gepasi  2 

Gepasi  2  is  a  free  software  package  developed  by  Pedro  Mendez  [23]  for  the 
intended  use  of  simulating  models  of  biological  and  chemical  reactions.  Its  de¬ 
velopment  was  based  on  biochemical  systems  theory  created  by  Savageau[29],  and 
metabolic  control  analysis  explained  by  Heinrich  et  al.  [15],  as  well  as  various  pre¬ 
existing  simulators.  It  has  an  easy  to  use  interface  and  will  run  on  MS-DOS  and 
MS-Windows.  It  has  no  plotting  capabilities,  but  has  the  capability  to  incorporate 
a  graphics  simulator.  It’s  purpose  is  to  mathematically  simulate  chemical  systems 
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with  both  predefined  and  user  defined  kinetic  equations.  It  conies  with  a  short  user 
manual  as  well  as  examples  of  simulations. 

Gepasi  2  first  requires  the  input  of  stoichiometric  eqiiations  representing  chem¬ 
ical  reactions.  It  has  the  capability  of  incorporating  up  to  45  metabolites  and  45 
reactions.  It  then  transforms  these  stoichiometric  equations  into  ordinary  differ¬ 
ential  equations  (ODE)  that  can  be  solved  using  a  ODE  solver.  Gepasi  2  uses  the 
free  solver  LSODA  developed  by  Petzold  [26].  LSODA  is  a  ordinary  differential 
equations  solver  that  can  test  whether  or  not  a  system  is  stiff  by  determining  the 
ratio  of  the  largest  and  smallest  eigenvalues  of  the  system.  If  the  ratio  is  large  in 
magnitude,  then  the  system  is  stiff  and  solved  using  a  stiff  solver.  If  the  ratio  is 
small  in  magnitude,  the  system  is  non-stiff,  a  non-stiff  solver  can  be  used.  LSODA 
solves  stiff  systems  using  the  backwards  differentiating  formula.  This  is  a  implicit 
method  whose  region  of  absolute  stability  contains  the  entire  negative  real  axis  and 
a  large  portion  of  the  complex  plane  so  as  to  retain  a  stable  solution.  LSODA 
solves  non-stiff  systems  using  the  Adams  method.  The  Adams  method  is  an  explicit 
method  with  a  smaller  region  of  stability  that  involves  time  steps.  It  requires  input 
from  the  user  to  define  its  time  step  size  by  asking  for  an  end  time  and  a  number  of 
intervals.  [2:pp  407-414] 

Other  functions  of  Gepasi  2  include  an  ability  to  solve  for  steady  states  and 
to  reduce  the  stoichiometric  matrix  defined  by  the  user.  It  uses  a  damped  Newton- 
Raphson  method  to  solve  for  steady  states  of  the  system.  As  a  precaution,  this 
method  ends  after  10^°  time  units  (in  seconds,  this  would  correspond  to  about  300 
years).  It  is  assumed  that  if  the  system  doesn’t  converge  at  that  point,  a  meaningful 
steady  state  cannot  be  reached.  It  reduces  the  stoichiometry  matrix  using  the  Gauss 
reduction  method. 

Gepasi  2  has  35  of  the  most  commonly  used  kinetic  equations  predefined  in 
the  simulator  so  as  to  be  easily  integrated  into  the  users  system.  These  kinetic 
types  include  the  Henri-Michaelis-Menton  equation  [34:ch  13]  and  the  Hill  equation 
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[34:pp  217-219]  for  both  forward  and  reverse  reactions.  It  also  has  the  capability 
of  implementing  kinetic  equations  designed  to  act  as  feedback  loops.  These  loops 
can  represent  either  activator  or  inhibitors  depending  on  the  kinetic  equation  used. 
User  defined  kinetic  types  can  also  be  implemented.  Mendes[23]  says  that  in  cases 
where  partial  derivatives  are  input,  they  can  be  solved  using  finite  differences. 

2.10  Other  Work 

2.10.1  V-Cell.  V-Cell  was  developed  by  Leslie  Loew  and  James  Schaff  of 
the  University  of  Connecticut  Health  Center.  Although  not  used  in  this  thesis,  it  is 
important  to  address  V-Cell  due  to  its  potential  as  a  simulator.  V-Cell  is  a  simiilator 
that  has  the  capability  to  model  the  cell’s  shape,  volume,  and  other  physical  features. 
It  also  models  how  molecules  diffuse  through  the  cell  by  using  a  system  developed 
through  the  understanding  of  actual  tests  of  diffusive  activity  in  cells.  Although 
V-Cell  is  a  smaller  simulator  than  E-Cell,  it  has  the  potential  for  creating  greater 
realism  within  a  cell. 

2.10.2  Global  Identifiably  Algorithms.  Once  cellular  functions  are  modeled 
and  formulated  as  first  order  ordinary  differential  equations,  they  need  to  be  solved. 
This  would  be  an  easily  solved  problem  if  the  rate  coefficients  were  known.  Unfortu¬ 
nately,  the  rate  coefficients  are  very  hard  to  measure.  Furthermore,  it  is  often  the 
case  that  small  errors  in  the  measurements  lead  to  large  errors  in  the  solution,  so 
the  measurements  must  be  extremely  accurate. 

Cobelli  et  al.[8]  considered  algorithms  such  as  the  Bunchberger  algorithm  to 
estimate  the  rate  constants  of  nonlinear  models.  They  incorporate  the  Bunchberger 
algorithm  into  a  larger  algorithm  to  check  for  global  identifiably,  a  property  that,  un¬ 
der  the  assumptions  of  noise-free  observations  and  error-free  model  structure,  guar¬ 
antees  uniqueness  of  the  solution. 
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2.11  Summary 

This  Chapter  described  several  examples  of  modeling  transcription  and  trans¬ 
lation  and  several  examples  of  modeling  metabolic  pathways.  Chapter  3  presents  a 
detailed  analysis  of  these  models  as  well  as  an  approach  for  analyzing  large  systems 
which  include  transcription,  translation,  and  metabolism  combined. 
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III.  Analysis  of  Models  of  Transcription,  Translation,  and 

Metabolism 

3. 1  Overview 

Models  describing  cellular  activity  have  been  created  to  predict  cellular  behav¬ 
ior  for  different  concentrations  of  component  molecules,  and  different  rate  constants. 
Many  of  these  models  are  described  using  systems  of  reaction  eqiiations  that  can  be 
reformulated  into  systems  of  differential  equations  that  track  contimious  movement 
of  each  component  though  time.  These  systems  of  differential  eqiiations  can  provide 
detailed  analysis  about  the  state  of  the  system.  This  chapter  will  begin  by  indicating 
how  to  linearize  about  a  steady  state  and  then  determining  local  behavior  from  the 
eigenvalues  of  the  linearized  system.  It  will  then  describe  models  of  equations  that 
exhibit  oscillatory  behavior  and  multiple  steady  states.  These  models  will  include  a 
linear  problem,  the  Brusselator  equation,  an  equation  developed  by  Schnackenberg, 
and  a  simple  multi-steady  state  problem.  This  chapter  will  conclude  by  discussing 
previously  developed  models  used  to  describe  transcription  and  translation.  These 
models  include  the  model  developed  by  Smolen  et  al.  and  incorporated  by  Hasty, 
Hasty’s  model,  and  Chens  larger  system  model.  By  better  understanding  these  mod¬ 
els,  it  will  be  easier  to  build  upon  and  reformulate  these  models  so  as  to  construct 
models  that  represent  Air  Force  needs. 

3.2  Analysis 

Mathematical  models  of  genetic  processes  can  evolve  into  systems  of  ordinary 
differential  equations.  When  evaluated,  these  systems  can  provide  detail  as  to 
how  the  model  behaves.  Techniques  for  analyzing  systems  of  ordinary  differential 
equations  have  been  available  for  many  years.  These  techniques  involve  finding  the 
steady  state  solution  and  checking  for  stability.  The  principal  reference  used  for 
this  analysis  is  Brauer  and  Nohel  [5].  They  lay  the  foundation  for  the  analysis  of 
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systems  of  autonomous  ordinary  differential  equations  by  giving  a  general  definition 
of  stability,  defining  stability  of  linear  systems,  and  explaining  how  to  analyze  the 
stability  of  nonlinear  systems  through  the  process  of  linearization.  Gray  and  Scott 
[11]  also  analyze  stability.  They  give  a  thorough  analysis  of  how  to  linearize  a 
complex  nonlinear  system  to  analyze  stability  inside  a  small  region.  This  analysis  is 
then  applied  to  the  metabolism  by  Heinrich,  Rapoport,  and  Rapoport,  as  an  example 
of  how  important  this  analysis  is  to  the  understanding  of  biological  systems. 

Brauer  and  Nohel  [5]  give  the  definition  of  both  stability,  asymptotic  stability, 
and  instability.  They  begin  by  defining  the  autonomous  system  below. 

y'  =  f(y)  where  y  is  a  n  dimensional  vector  (3.1) 

They  then  define  the  steady  state  point  which  they  refer  to  as  the  critical  point. 

Let  y  =  yo  be  a  critical  point.  Then  f(yo)  =  0  (3.2) 

If  there  exist  no  other  critical  points  within  a  neighborhood  of  yo,  then  yo  is  an 
isolated  critical  point. 

Next,  they  outline  the  definition  of  stability  followed  by  the  definition  of  as¬ 
ymptotic  stability. 

Definition  1  The  steady  state  solution  yo  of  equation  (3.1)  is  said  to  he  stable  if  for 
eaeh  number  £  >  0  we  ean  find  a  number  6{e)  >0  such  that  ifipit)  is  any  solution  of 
equation  (3.1)  hawing  ||'0(to)  ^  yo||  <  b,  then  the  solution  'ip{t)  exists  for  all  t  >to. 

Definition  2  The  steady  state  solution  yo  is  said  to  be  asymptotically  staMe  if  it 
is  staMe  and  if  there  exists  a  number  5(e)  >  0  such  that  if  1^(1)  is  any  solution  of 
equation  (3.1)  having  ||'0(to)  ^  yo||  <  b,  then  limt^oo'0(t)  =  yo. 

Definition  3  The  steady  state  solution  yo  is  said  to  he  unstaMe  if  it  is  not  stuble. 
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Brauer  and  Nohel  then  explain  that  the  simplest  general  systems  for  which  sta¬ 
bility  is  completely  defined  is  the  linear  system.  Given  the  linear  system  in  equation 
(3.3)  below,  were  A  is  a  real  constant  n  x  n  matrix,  stability  can  be  determined  by 
considering  the  eigenvalues  of  A. 


y'  =  Ay 


(3.3) 


The  specific  criteria  for  stability  about  the  origin  is  given  in  the  following 
theorem. 


Theorem  1  “If  all  eigenvalues  of  A  have  nonpositive  real  parts  a,nd  all  those  eigen¬ 
values  with  zero  real  parts  are  simple,  then  the  solution  y  =  0  of  (3.3)  is  staMe.  If 
and  only  if  all  eigenvalues  of  A  have  negative  real  parts,  the  zero  solution  of  (3.3)  is 
asymptotically  stable.  In  fact,  in  this  case  if'ip(t,tQ)  denotes  the  fundamental  matrix 
of  (3.3)  which  is  the  identity  a,t  t  =  tQ,  then  'ip{t,  to)  =  exp{{t  —  to) A)  and  there  exist 
constants  K  >  0,  a  >  0  such  that  |■0(^,  to)|  <  K e,xp{—cr{t  —  to))  where  {to  <t  <  oo) 
with  (7  >  0  in  the  case  thal  all  eigenvalues  of  A  have  negalive  real  parts  and  a  =  0 
if  there  are  simple  eigenvalues  with  zero  real  part.  If  one  or  more  eigenvalues  of  A 
hare  a  positive  real  part,  the  zero  solution  of  (3.3)  is  unstable.”  [5:p  151] 

This  analysis  of  linear  systems  is  useful  in  describing  localized  steady  states 
stability  of  some  nonlinear  systems.  This  is  done  though  the  process  of  linearization. 
Gray  and  Scott  [11]  followed  the  definition  provided  by  Brauer  and  Nohel  [5]  by 
explaining  that  local  asymptotic  stability  is  defined  by  the  behavior  of  a  system 
very  near  the  steady  state  point.  If  a  small  perturbation  from  steady  state  grows, 
the  steady  state  is  unstable.  Where  as  if  it  decays,  the  system  is  stable.  The 
qualification  ‘local’  implies  the  system  may  exhibit  different  stability  behavior  for 
larger  perturbations.  They  then  explain  that  linearization  is  a  perturbation  method 
that  allows  nonlinear  systems  to  be  treated  locally  as  linear  systems.  Brauer  and 
Nohel  [5]  give  a  general  description  of  how  to  linearize  a  two  dimensional  system 
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about  a  steady  state  point  located  at  the  origin.  Gray  and  Scott  give  a  more  explicit 
derivation  of  how  to  linearize  the  two  dimensional  system  about  any  steady  state 
point.  This  description  is  expanded  upon  by  Heinrich,  Rapoport,  and  Rapoport  [15] 
into  an  n  dimensional  case.  However,  Heinrich,  Rapaport,  and  Rapaport  give  little 
detail  leading  to  their  results.  The  detail  of  the  n  dimensional  case,  is  given  below. 

First,  consider  the  nonlinear  system  given  by  eqiiation  (3.1).  Assume  the 
system  has  a  steady  state  value  at  f(yss)  =  0.  Now,  a  small  perturbation  away 
from  the  steady  state  is  given  below. 


y  =  Yss  +  Ay  where  \Ayi\  <<  1  for  (i  =  l...n)  (3.4) 

Then  the  rate  of  change  of  the  perturbation  Ay  is  obtained  by  substituting 
equation  (3.4)  into  equation  (3.1).  Note  that  y^^  is  a  constant  so  =  0. 


d(y,,  +  Ay)  djyj  ^  d{Ay) 

(it  dt  dt 


d{Ay) 

dt 


f(y.!  +  Ay) 


(3.5) 


Using  a  Taylor  Series  Expansion  about  the  steady  state  point,  this  equation 
expands  to  the  following. 


=  /i(yss)  +  'Yl  terms  for  (j  =  l...n)  (3.6) 

Now  fjijss)  =  0  for  all  j  and  (higher  order  terms  =  0)  by  perturbation  theory 
coupled  with  the  assumption  |Aj/j|  <<  1  for  {i  =  l...n)  so  equation  (3.6)  simplifies 
to  the  following. 


dt 


E 


dfjiyss) 

% 


A?/i 


for  (j 


(3.7) 
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Now  by  introducing  the  Jacobian  matrix  J,  the  above  equation  is  rewritten  in 
vector  form  below. 


^(Ay) 

(it 


where  J(y) 


J{yss)Ay 

for  (i  = 

(  S/i(y) 

dfi  (y) 

9  My) 

\ 

dy\ 

dy2 

9yn 

9/2  (y) 

9/2  (y) 

; 

dyi 

9j/2 

dfniy) 

dfn(y) 

) 

\  dyi 

dyn 

(3.8) 


(3.9) 


Furthermore,  Jiyss)  is  a  constant  matrix,  so  equation  (3.8)  is  similar  to  equa¬ 
tion  (3.3),  and  theorem  1  applies.  The  eigenvalues  of  the  system  above  are  denoted 
as  A  in  the  system  below. 


Det  I  A/  —  J(yss)|  =  0  where  /  is  an  n  x  n  Identity  matrix  (3.10) 

Solving  for  A  and  applying  theorem  1  defines  the  local  steady  state  stability  for  the 
nonlinear  system. 

Gray  and  Scott  [11]  focused  on  the  linearization  of  a  two  dimensional  model 
to  evaluate  various  types  of  stability  while  retaining  a  mathematical  base.  The 
linearized  version  of  their  two  dimensional  system  is  as  follows. 


d{Aa)  dfss  .  dfss  .  ^  ^ 

=  ^Aa+^A/l  =  aiiAa  +  ai2A/l 
(it  da  op 

(i(A/3)  dgss  ^  dgss  ^  ^  a  i  \  q 

— 1—  =  Aa-f-^— A/l  =  a2iAa  +  a22A/l 

(it  da  dp 


(3.11) 


Then,  using  equation  (3.10)  to  solve  for  the  eigenvalues,  the  quadratic  equation  below 
is  formed 
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—  tr{J)X  +  det(  J)  =  0 


(3.12) 


where  J  =  \  1 ,  tr{J)  =  an  +  022,  and  det( J)  =  011022  —  O12O21. 

\  O21  O22  J 

Now,  equation  (3.12)  has  the  solutions  Ai,2  =  |{i?’(d^)  ±  -ytr(Jp'^^Tdet(J)}. 
The  different  possible  eigenvalue  solutions  provide  for  a  mimber  of  different  types  of 
stability  scenarios.  These  scenarios  along  with  the  requirements  for  their  existence 
are  listed  in  the  book  Chemical  Oscillations  and  InstaMlities  by  Gray  and  Scott 
[ll:pp  65-68]. 

3. 3  Specific  cases  of  staMlity 

In  this  section  simple  mathematical  models  that  exhibit  different  characteristics 
of  stability  are  considered.  The  first  model  is  a  simple  two  component  system  that 
displays  basic  characteristics  of  stability  that  are  sensitive  to  different  parameter 
values.  A  sensitivity  analysis  is  performed  on  a  parameter  within  the  model  to 
examine  its  effect  on  the  models  stability.  Also  considered  are  two  models  of  reaction 
equations  that  exhibit  oscillatory  steady  state  behavior.  These  models  include  the 
Brusselator  equation,  and  an  equation  developed  by  Schnackenberg.  The  last  model 
included  will  consist  of  a  small  system  that  shows  interesting  behavior  by  presenting 
multiple  steady  states  with  different  parameter  values.  By  looking  at  some  basic 
mathematical  models  that  define  different  types  of  stability  and  then  incorporating 
them  into  biological  systems  to  achieve  similar  stability  characteristics,  more  can  be 
understood  about  the  biological  system  as  a  whole. 

3.3.1  A  two  component  system.  Consider  the  following  two  component 
system. 


fcs 


X  4 

fci 


Ike 

Y 


(3.13) 
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This  is  a  small  system  that  can  be  involved  in  any  larger  biological  system.  It 
can  be  represented  as  a  set  of  differential  equations  in  the  following  manner. 


(it 

m 

(it 


h[Y]  -  UX]  -  h[X]  +  h 
hlX]  -  HY]  -  k^[Y\  +  h 


(3.14) 

(3.16) 


This  set  of  differential  equations  can  be  simplified  to  the  following  set  of  eqiia- 

tions. 


i\X\ 

(it 

m 

(it 


h[Y]  -  k;[x]  +  k, 
k2[X]  -  kl[Y]  +  A)6 


(3.16) 


where  k^  =  k2  +  k^  and  kl  =  ki  +  k^ 

Now,  the  stability  associated  with  the  above  equations  can  change  dramatically 
given  small  changes  to  specific  rate  constants.  For  example,  say  the  rate  constants 
are  defined  as  follows,  where  there  is  a  small  perturbation  that  changes  the  rate 
constants  in  a  manner  defined  by  e. 

With  k,  =  3,k2  =  ^3*  =  3  +  £,  kl  =  l,ks  =  0,ke  =  0,x=  [X],  y  =  [T], 


then  system  (3.16)  becomes 


and  J  = 


-(3  +  e)  3 

(2+g)" 


with  a  steady  state  point  at  x 


dy_ 

dt 

0  and  y 


JV 

0. 


(3.17) 


(3.18) 
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The  solution  for  this  system  takes  the  general  form 

X  =  (3.19) 

Y  =  cae^it  +  C4e^2t  ^^20) 

where  Ai  and  Aa  are  the  eigenvalues  of  equation  (3.18),  and  are  determined  using 
equation  (3.12)  to  be  Ai,a  =  ^{tr(J)  ±  ^tr(J)^  -  4det(J)},  where  tr(J)  =  -(4  +  e), 
and  det( J)  =  (3  +  e)  —  3^^^^  =  — 2e  — 

Thus,  Ai,2  =  ^{-(4  +  e)  ±  (2  +  e)'}  =  -^(4  +  e)  ±  (2  +  e)  (3.21) 

So  Ai  =  -e,  and  Aa  =  — 4  —  -e  (3.22) 

2  2 

In  this  system,  e  is  restricted  to  perturbations  no  less  than  —3  because  the 
kinetic  constant  is  restricted  to  positive  values.  Therefore,  within  the  interval 
that  e  is  restricted  to,  there  exists  two  bifurcation  points.  The  first  is  at  e  =  0,  and 
the  second  is  at  £  =  —  |.  At  each  bifurcation  point,  the  system  changes  stability. 
When  £  is  within  the  interval  [— 3,  — |),  Ai  <  0,  and  Aa  >  0.  So,  the  steady  state 
point  acts  as  an  unstable  saddle  point.  When  e  is  within  the  interval  (— 1,0), 
Ai  <  0,  and  Aa  <  0.  So,  the  system  converges  asymptotically  to  the  steady  state 
point.  When  £  is  within  the  interval  (0,  oo),  Ai  >  0,  and  Aa  <  0.  So  again,  the 
steady  state  point  acts  as  an  unstable  saddle  point.  This  example  shows  that  a 
small  change  in  a  single  parameter  can  change  the  stability  of  a  system  greatly.  In 
cases  where  parameters  are  a  bit  relaxed,  or  changed  due  to  an  outside  influence,  a 
stability  analysis  is  an  important  tool  in  describing  the  system. 
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3.3.2  The  Brusselator.  Limit  cycle  behavior  can  exist  in  both  chemical 
systems  [9:pp  2352-2357]  and  biological  systems  [15].  For  this  reason,  it  is  necessary 
to  consider  the  modeling  techniques  associated  with  limit  cycle  oscillatory  behavior. 
First  proposed  by  Prigogine  and  Lefever  [27]  in  1967,  the  Brusselator  is  a  popular 
chemical  model  [9]  that  describes  limit  cycle  oscillatory  behavior.  It  is  given  by  the 
following  chemical  equations. 


A  - 

X 

(3.23) 

B-hX  - 

Y  +  D 

(3.24) 

2X-hr  - 

3X 

(3.25) 

X  - 

E 

(3.26) 

These  chemical  equations  can  be  transformed  into  a  system  of  differential  equa¬ 
tions  using  the  techniques  summarized  in  Chapter  2  and  presented  by  Heinrich, 
Rapoport,  and  Rapoport  [15].  First,  consider  the  following  rate  equations  describ¬ 


ing  the  reaction  equations  (3.23)-(3.26). 

Vi  = 
V2  = 
Vs  = 

V4  = 


/i| 

(3.27) 

(3.28) 

(3.29) 

:^i 

(3.30) 

These  can  be  used  to  describe  the  rate  of  change  of  each  component  as  it  relates  to 
the  entire  model.  This  results  in  the  following  set  of  differential  equations. 
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-Vi  +  ^2  +  2^3  - 


(3.31) 

(3.32) 


d[X] 

dt 


m 

dt 


=  -V2  +  V3 


By  assuming  [A]  =  a,  [B]  =  b,  [X]  =  x,  and  [y] 
this  set  of  equations  simplifies  to  the  following. 


y  where  a  and  b  are  held  constant, 


dx 

dt 

dy 

dt 


a  —  {b  +  l)x  +  x^y 
bx  —  x?y 


(3.33) 

(3.34) 


The  system  is  in  steady  state  when  ^  =  0  and  ^  =  0.  Therefore,  solving  for  the 
steady  state  solution  involves  simply  solving  the  following  system  for  x  and  y. 


0  =  a  —  (h-\- l)x  +  x'^y  (3.35) 

0  =  bx  —  x^y  (3.36) 

From  equation  (3.36),  y  =  f  •  Substituting  y  into  equation  (3.35)  yields  x  =  a, 
which  yields  y  =  ^-  Therefore,  the  steady  state  point  is  as  follows. 

X  =  a  (3.37) 

y  =  I  (3-38) 

As  a  result  of  linearizing  equations  (3.33)  and  (3.34),  the  local  stability  around 
the  steady  state  can  be  determined.  From  equation  (3.11), 
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J  = 


2xy  —  (6  +  1) 
b  —  2xy  — 


Jss 


b-1 


The  eigenvalues  are  calculated  by  solving  the  characteristic  equation  3.39. 


Det  |AJ  —  Jss\  =  Det 

Thus,  the  eigenvalues  are  A 

6  =  +  1, 


A 


1,2 


b+1 

b 


A  + 

-l-a2 


=  A^  +  A(a2  +  1  -  &)  +  (3.39) 


±  v^(o2  +  1  -  &)2  -  4a2}.  By  setting 


Ai^2  =  i\/— (2^  =  zizia,  (3.40) 

and  Ai,2  are  both  purely  imaginary.  This  is  known  as  a  Poincare-Andronov-Hopf 
bifurcation  point  [12]  because  any  perturbation  in  h  such  that  +  1  ±  e  for 

e  >  0  will  lead  to  either  stability  or  divergent  instability.  This  can  be  tested 
through  stability  analysis.  If  6  =  +  1  —  then  Ai^2  =  ±  \Je'^  —  da^}  which 

leads  to  a  damped  oscillatory  approach  to  steady  state  due  to  the  negative  real  part 
of  A.  If  6  =  +  1  +  £,  then  Ai,2  =  |{£  ±  \/e^  —  4a^}  which  leads  to  oscillatory 

divergence  due  to  the  positive  real  part  of  A. 

3.3.3  An  Example  developed  by  Schnackenberg.  Schnackenberg[30]  sys¬ 
tems  are  simple  chemical  oscillators.  Consider  the  following  model  presented  by 
Schnackenberg. 
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2X  +  F  3X 

(3.41) 

A^Y 

(3.42) 

X  A  B 

k-3 

(3.43) 

The  velocity  though  each  chemical  equation  is  given  as, 


=  -HXf[Y]  (3.44) 

V2  =  -k2[A] 

U3  =  -h[X]  +  k_^[B] 

The  differential  equations  corresponding  to  (3.41)- (3.43)  are 

2vi  -  3ui  +  U3  =  -vi  +  ^3  (3.45) 

Vi  -  V2 

Using  equation  (3.44)  and  letting  [X]  =  x,  [T]  =  y,  [A]  =  a,  and  [B]  =  /3 
yields  the  following. 


(it 

m 

(it 


dx 

dt 

% 

dt 


kix^y  —  k^x  +  k^sP 


—k\x^y  +  k20t 


(3.46) 
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Nondimensionalizing  equation  (3.46)  by  setting  x  =  y  =  Yq-i],  and  t  =  Tqt 
yields  the  following. 


dr 

dr] 


hXoYoT^ev  -  hToC  + 

7  v2rr  t2  ,  k2Toa 
-hiX^To^  r]+—— 


(3.47) 

(3.48) 


IfTn  = 


"^0  “  \i  ^  ~  ^0’  then  kiXoYQTo  —  1,  hiX^T^  —  1,  k^To  —  1 


0 

and  equations  (3.47)  and  (3.48)  become 


dr 

dr] 

dr 


ev-^+b 

-i^r]  +  a 


(3.49) 

(3.50) 


where  b  =  ^  and  a  = 

The  steady  state  can  be  found  by  setting  ^  =  0,  and  ^  =  0,  and  then  solving 
the  following  system  for  ^  and  rj  in  terms  of  a  and  h. 


0  =  ^^V-^  +  b  (3.51) 

0  =  -^2^  +  a  (3.52) 

By  adding  equation  (3.52)  to  equation  (3.51),  the  equation  0  =  —^  +  b  +  a,  can 
be  solved  in  terms  of  ^  to  obtain  ^  =  o  +  6.  Then  substituting  ^  back  into  equation 
(3.52),  yields  0  =  —  (a  +  6)^r/  +  a,  which  yields  rj  =  Thus,  the  steady  state 

solution  is 
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(3.53) 

(3.54) 


^  =  a  +  b 
^  (a  +  bf 

Then,  linearizing  equations  (3.49)  and  (3.50),  using  eqiiations  (3.11)  yields. 


J 


e  \  ^ 

-2^V  ^ 


2a 

(a+b) 


—2a 

(a+b) 


{a  +  bf  \ 
-{a  +  bf  ) 


(3.55) 


Solving  the  characteristic  equation, 


Det  I XI  —  Jss 


A 


Det 


2a 

(a+6) 


+  1 


2a 

{a+b) 


+  A(((i  +  b)^  +  1 


-(a +  5)2 
A  +  (u  +  5)2 


(3.56) 


yields  the  eigenvalues  Ai, 2  =  5{-(a+5)2-l+^^±y^((a  +  5)^  +  1  -  -  4(a  +  5)2} 

By  setting  (a  +  5)2  +  1  —  =  0,  the  system  has  a  Poincare- Andronov-Hopf 

bifurcation  point  where  (a  +  5)^  =  a  —  5  because  the  real  part  of  the  eigenvalues 
disappears  thereby  leaving  purely  imaginary  eigenvalues  and  therefore,  a  limit  cycle. 

For  a  small  perturbation  from  the  bifurcation  point  in  one  direction,  the  steady  state 
acts  as  a  spiral  sink  due  to  a  small  negative  real  part  added  to  the  eigenvalues,  while 
for  a  perturbation  in  the  other  direction,  the  bifurcation  point  acts  as  a  spiral  source 
due  to  a  small  positive  real  part  added  to  the  eigenvalues. 

In  evaluating  the  simple  Schnackenberg  model,  there  exist  many  similarities 
to  the  Brusselator  model.  Both  models  include  an  autocatalytic  chemical  equation. 
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both  models  measure  the  rate  of  change  of  two  species  in  very  small  systems,  and 
both  models  have  limit  cycles.  If  the  Brusselator  were  mathematically  equivalent 
to  Schnackenbergs  model,  then  there  would  be  a  one-to-one  onto  mapping  from  the 
Brusselator  model  to  Schnackenbergs  model  defining  an  isomorphism.  However, 
the  models  are  distinctly  different.  Schnackenberg  states  that  his  oscillator  is  “not 
isomorphic  to  the  Brusselator  and  chemically  originates  from  a  different  and  smaller 
system  of  reactions.”  [30] 

3.3.4-  ^  small  system  involving  bistability.  Under  specific  conditions,  a 

system  can  have  multiple  steady  states.  This  sometimes  translates  to  miiltiple  stable 
steady  states.  The  next  system  “shows  multiple  steady  states  and  has  been  proposed 
by  Schlogl  (1972)  as  a  particularly  simple  and  analytically  solvable  example.”  [30] 
It  is  interesting  because  its  simplicity  allows  it  to  be  involved  in  larger  systems.  The 
system  is  as  follows. 


2X  +  A  ^3X 

(3.57) 

k-i 

k2 

X  A  B 

k-2 

(3.58) 

The  velocity  of  a  single  component  through  these  chemical  equations  is  given 

below 


=  -ki[X]^[A]  +  k_i[Xf  (3.59) 

V2  =  -k2[X]  +  k-2[B] 

This  allows  for  each  component  to  be  tracked  through  the  model  allowing  a 
system  of  differential  equations  to  be  created. 
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(It 


2vi  -  3vi  +V2  =  -Vi  +  V2 


(3.60) 


Now,  substituting  equation  (3.59)  into  equation  (3.60),  the  rate  of  change  of 
the  one  component  system  can  be  tracked  explicitly. 


d[X] 

(It 


=  h[Xf[A]  -  k.r[Xf  -  k2[X]  +  k.2[B] 


(3.61) 


To  nondimensionalize,  let  [X]  =  Xqx,  t  =  Tqt,  [A]  =  A,  and  [B]  =  B  so  that 


^  =  kiXoTo[Ay  -  k.iXlTox^  -  k2Tox  + 

(IT  .Aq 

If  To  =  -^  and  Xq  =  J then  k^iX^Tf)  =  1,  k2TQ  =  1  and 


(3.62) 


dx  2  s 

—  =  ax  —  X  —  X  -\-  0 

dr 


(3.63) 


where  a  =  kiXoTo[A],  and  b  = 

Solving  for  steady  state,  leads  to  a  cubic  polynomial  with  three  roots. 


x^  —  ax?  +  X  —  6  =  0 


(3.64) 


Understanding  oscillatory  behavior  and  multiple  steady  state  behavior  in  math¬ 
ematical  systems  may  lead  to  an  understanding  of  how  to  detect  and  construct  the 
behavior  in  biological  models.  With  this  knowledge,  it  may  be  possible  to  gain 
more  understanding  of  a  large  biological  systems  by  incorporating  these  models  that 
represent  small  systems  into  them. 
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3-4  Smolen  et  al.  and  Hasty  et  al. 


Smolen,  Baxter,  and  Byrne[32]  developed  a  model  to  test  the  dynamic  prop¬ 
erties  of  a  genetic  regulatory  system.  They  proposed  a  model  that  ’’captures  the 
salient  features  of  transcription  factors,  dimerization,  binding,  and  phosphorylation- 
dependent  regulation  of  transcription.”  They  evaluated  the  system  to  find  the 
steady  state  concentration  of  the  transcription  factor.  They  then  plotted  the  steady 
state  concentration  of  the  Transcription  Factor  with  respect  to  a  parameter  a  to 
show  how  changes  in  a  generate  bistability  in  the  system.  Unfortunately,  Smolen, 
Baxter,  and  Byrne  neglect  to  give  any  sort  of  derivation  to  show  how  their  model 
leads  to  their  conclusions.  Hasty  et  al.[14]  picks  up  on  the  conclusions  of  Smolen, 
Baxter,  and  Byrne  by  using  the  model  to  evaluate  additive  noise.  Hasty  et  al.  begins 
by  symbolically  writing  out  the  model  in  the  form  of  a  system  of  reaction  equations. 
However,  from  there.  Hasty  et  al.  jump  directly  to  form  their  conclusions  without 
including  any  derivation  in  their  model.  It  is  important  to  understand  the  deriva¬ 
tion  of  the  model  in  order  to  judge  its  validity.  Furthermore,  understanding  how 
to  construct  bistability  in  this  model  will  help  to  formulate  models  that  represent 
future  Air  Force  needs. 

The  construction  of  this  model  begins  with  Hasty  et  al.’s  system  of  reaction 
equations. 


ki 

2X  ^  X2 

(3.65) 

k-i 

D  +  X2^  DX2 

(3.66) 

k—2 

DX2  +  P^  DX2  +  nX 

(3.67) 

xhA 

(3.68) 
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This  system  consists  of  two  fast  reactions,  (3.65)  and  (3.66),  that  represent 
the  binding  of  two  X  proteins  to  form  a  dimer  molecule  that  then  binds  to  a  DNA 
promoter  site  to  form  a  DNA  promoter  complex.  Then,  with  the  introduction  of 
RNA  polymerase,  the  third  equation,  represents  both  transcription  and  translation 
that  takes  place  at  rate  kf  to  create  n  X  proteins.  Finally,  the  fourth  equation 
represents  degradation  of  protein  X  at  the  slow  rate  of  kd.  Together,  these  equa¬ 
tions  represent  a  simple  model  of  protein  regulation  through  gene  expression.  An 
additional  equation  not  listed  by  Hasty  et  al.  must  also  be  included  to  reproduce 
the  results  presented  in  their  paper.  Equation  (3.69)  represents  “the  basal  rate  of 
production  of  protein  X,  i.e.,  the  low  baseline  expression  rate  in  the  absence  of  a 
transcription  factor.”  [14]  It  will  be  determined  if  this  basal  rate  of  prodiiction  is 
necessary  for  the  bistability  in  Chapter  4. 

INPUT  ^  X  (3.69) 

These  equations  can  be  transformed  into  a  system  of  differential  equations  using 
the  techniques  outlined  by  Heinrich,  Rapoport,  and  Rapoport  [15]  and  displayed  in 
chapter  2.  In  the  equations  below,  Uj  represents  the  flux  of  a  single  substrate 
molecule  through  one  of  the  chemical  equation  described  above.  For  example,  in 
equation  (3.70),  Vi  represents  the  reversible  flow  of  one  molecule  from  substrates  to 
products  at  a  rate  of  ki  and  back  again  with  a  rate  of  k-i. 


■^1 

=  k^^{x^]-h[xf 

(3.70) 

V2 

=  k^2[DX2]-k2[D\[X2] 

(3.71) 

V3 

=  -k3[DX2][P\ 

(3.72) 

V4 

=  -k^[X] 

(3.73) 

Vh 

=  r 

(3.74) 
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These  equations  can  then  be  used  to  track  the  flow  of  each  molecule  through 
the  whole  system.  Again,  this  was  described  by  Heinrich,  Rapoport,  and  Rapoport 
[15]  in  Chapter  2.  Together,  the  production,  reaction,  and  degradation  of  each 
molecule  is  tracked  through  the  reaction  equations  using  eqiiations  (3.70)  through 
(3.74)  and  described  with  the  differential  equations  presented  next.  For  example, 
two  X  substrate  molecules  can  be  tracked  through  eqiiation  (3.65).  They  are  then 
lost  through  reaction,  but  recreated  n  times  over  via  eqiiation  (3.67),  and  recycled 
though  the  system  until  they  degrade  in  equation  (3.68).  At  the  same  time,  X 
molecules  are  being  created  at  a  basal  rate  described  by  eqiiation  (3.69). 


i\X] 

(it 

4X2] 

(it 

dlD] 

(it 

(i[DX2] 

(it 


2vi  -  nvs  +  U4  +  U5 


-Vi  +  V2 


V2 


-V2  +  U3  -  U3 


(3.75) 

(3.76) 

(3.77) 

(3.78) 


The  degraded  state  of  protein  X  that  takes  the  form  A,  is  not  an  active  part 
of  the  system,  and  therefore,  not  considered  as  a  rate  that  drives  the  system.  Fur¬ 
thermore,  the  concentration  of  RNA  polymerase  P  is  considered  constant  over  time, 
[P]  =  po,  so  its  rate  of  change  is  not  included  in  the  system  above.  Equations  (3.75) 
through  (3.78)  are  then  expanded  to  yield  the  rate  equations. 
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d[X] 

dt 

d{X2] 

dt 

d{D] 

dt 

d{DX^] 

dt 


2(fc_i[X2]  -  h[X]^)  +  nh[DX2]po  -  h[X]  +  r 
-ik_^[X2]  -  h[Xf)  +  k_2[DX2]  -  k2[D][X2] 
k_2[DX2]  -  k2[D][X2] 

-{k.2[DX2]  -  k2[D][X2]) 


(3.79) 

(3.80) 

(3.81) 

(3.82) 


Now,  (3.81)  and  (3.82)  are  negatives  of  each  other.  This  leads  to 


d[D]  ^  d[DX2]  _  d{[D]  +  [DX2]) 
dt  dt  dt 

or  more  appropriately 


(3.83) 


[D]  +  [DX2]  =  dr  where  dr  is  a  constant  (3.84) 

Equation  (3.84)  says  the  amount  of  free  DNA  promoter  sites,  [D],  and  the 
amount  of  DNA  promoter  complexes,  [DX2],  makes  up  the  total  number  of  DNA 
promoter  sites,  dr,  in  the  system. 

At  this  point,  Hasty  et  al.’s  [14]  system  of  reaction  equations  has  been  trans¬ 
formed  into  a  system  of  differential  equations  that  represent  the  model  created  by 
Smolen,  Baxter,  and  Byrne[32].  From  here,  using  appropriate  assumptions  about 
the  reaction  speeds  of  the  chemical  equations,  (3.65)  though  (3.68),  Hasty  et  al. 
were  able  to  simplify  the  equations,  (3.79)  through  (3.84),  down  to  only  one  equa¬ 
tion  who’s  rate  was  determined  only  by  parameter  constants  and  concentrations  of 
protein  X.  The  assumption  was  that  equations  (3.65)  and  (3.66)  “are  orders  of 
magnitude  faster  than  transcription  and  degradation.”  [14]  The  assumption  that 
equations  (3.65)  and  (3.66)  are  fast  reactions  implies  that  the  contents  whose  concen¬ 
trations  change  only  within  the  confines  of  those  equations  can  be  treated  as  steady 
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state.  This  assumption  is  important  because  the  concentrations  of  both  X2  and 
D  are  represented  only  in  equations  (3.65)  and  (3.66).  Furthermore,  although  the 
concentration  of  DX2  is  represented  in  the  slow  eqiiation  (3.67),  its  concentration 
remains  constant  through  that  equation  and  only  fluctuates  with  eqiiation  (3.66). 
Thus,  =  0,  =  0,  and  =  0  because  the  concentrations  of  X2,  H,  and 

DX2  are  all  assumed  to  exhibit  steady  state.  With  this  new  information,  equations 
(3.79)  through  (3.84)  become 


^  =  2{k_^[X2]-h[Xf)  +  nks[DX2]po-h{X]+r  (3.85) 

0  =  -{k.i[X2]  -  ki[Xf)  +  k.2[DX2]  -  k2[D][X2]  (3.86) 

0  =  k.2[DX2]  -  k2[D][X2]  (3.87) 

(It  =  [D]  +  [DX2]  (3.88) 


This  system  can  be  simplified  in  the  following  manner.  Solving  equation  (3.87) 
for  DX2  in  terms  of  D  and  X2  yields 

[DX2\  =  ^{D\{X2\  (3.89) 

Then  substituting  equation  (3.89)  into  equation  (3.86)  results  in 

0  =  -(A;_i[X2]  -  k^[Xf)  +  k_2^[D\[X2]  -  k2[D][X2\  (3.90) 

k-2 

Solving  equation  (3.90)  in  terms  of  X2  yields 


IV2I  =  X[xp  (3.91) 

Substituting  equations  (3.89)  and  (3.91)  into  equation  (3.85)  and  simplifying  results 
in 
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(3.92) 


^  =  ’'P«-r-^k3lD][Xf  -  h[x\  +  r 

(it  fC_i  K—2 

Equation  (3.88)  can  then  be  used  to  remove  the  dependence  on  [D],  First  substitute 
equation  (3.89)  followed  by  equation  (3.91)  into  equation  (3.88)  to  obtain 

dr  =  [0|  +  (3.93) 

Solving  for  [D]  yields 


[D]  = 


1  +  ifcfcW' 


(3.94) 


This  equation  can  then  be  inserted  into  equation  (3.92)  to  eliminate  the  de¬ 
pendence  on  [D]. 


d[X\  np„-t-,f^,hdTlXp  _ 

dt  ■  l  +  ^ 


(3.95) 


Finally,  by  defining  =  ^,  a  =  np^KJ<2hdT,  P  =  K^K2, 

'j  =  k4,  6  —  r,  and  u  =  [X],  the  rate  equation  for  the  protein  X  can  be  displayed  in 
the  form  which  Hasty  et  al.[14]  and  Smolen  et  al.  [32]  present  it. 


du 

dt 


au^ 

1  +  PvP- 


—  ^u-\-  8 


(3.96) 


Hasty  et  ah,  and  Smolen  et  al.  then  examined  how  the  steady  state  value  of  u 
changed  as  a  function  of  a.  They  incorporated  a  bifurcation  plot  to  track  the  steady 
state  concentration  of  u  with  respect  to  a,  and  discovered  that  as  a  increases  over 
a  particular  interval,  u  attains  multiple  steady  states  and  in  fact  becomes  bistable. 
Again,  neither  Hasty  et  al.  or  Smolen  et  al.  give  explicit  derivation  showing  how  they 
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achieved  their  results.  Therefore,  we  will  now  analyze  steady  state  concentration  of 
u  with  respect  to  a  in  an  attempt  to  recreate  their  results. 

First,  set  %  =  0  in  order  to  evaluate  the  steady  state  concentration  of  u. 


0 


av^ 

1  +  jSv? 


—  yn  +  (5 


(3.97) 


Then,  multiply  through  by  the  denominator  and  simplify  to  obtain  the  following 
cubic  polynomial. 


0  =  n^ 


((5/3  +  tt)  2  ,  7  ^ 

R  “  a 

IP  IP  IP 


(3.98) 


Then,  for  different  parameter  settings,  this  equation  can  be  solved  for  u.  For 
example,  Hasty  set  /3  =  1,  7  =  10,  and  6  =  A  and  then  solved  for  u  with  respect 
to  ot  over  the  interval  a  =  [0,70].  With  this  choice  of  parameters,  equation  (3.98) 
becomes 


0  =  ^  _  04  (3.99) 

For  each  value  of  a,  the  roots  of  this  polynomial  can  be  determined. 

3.5  Hasty  et  a,l. 

The  Air  Force  is  searching  for  rapid  ways  to  construct  models  of  cellular  sys¬ 
tems.  They  are  interested  in  understanding  the  unique  genetic  responses  that  come 
from  novel  Air  Force  chemicals  introduced  into  cellular  systems.  One  way  to  quickly 
create  a  model  that  incorporates  desired  responses  is  by  considering  other  related 
models.  Hasty  et  al.[13]  shows  that  it  is  possible  to  take  a  genetic  model  that  incor¬ 
porates  normal  cellular  activity  and,  using  it  as  a  template,  construct  a  new  model 
that  incorporates  an  elaborate  cellular  response  to  outside  stimuli.  Understanding 
how  Hasty  et  al.  drew  from  the  previous  model  developed  by  Smolen  [32]  to  create  a 
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new  model  will  help  the  Air  Force  to  quickly  generate  it’s  own  models  using  similar 
techniques. 

Hasty  drew  from  the  model  constructed  in  the  previous  section  to  create  a 
model  representing  two  different  tracks  the  bacteriophage  A  virus  may  take  to  infect 
a  cell.  Hasty  et  al.  began  by  writing  out  the  reaction  eqiiations  that  represent  the 
cellular  processes  taking  place.  One  of  the  assumptions  in  [13]  is  that  eqiiations 
(3.101)  though  (3.104)  represent  fast  reactions,  where  as  the  others  represent  slow 
reactions. 


Input  A  X 

(3.100) 

ki 

2X  ^  X2 

k-i 

(3.101) 

D  +  X2^  DX2 

k-2 

(3.102) 

D  +  X2^  dx* 

k-3 

(3.103) 

DX2  +  X2^  DX2X2 

k—4 

(3.104) 

DX2  +  P^  DX2  +  P  +  nX 

(3.105) 

X 

(3.106) 

Equation  (3.100)  is  not  specifically  included  in  Hasty  et  al.’s  system  of  chemical 
equations,  but  does  represent  the  basal  rate  of  production  of  X  as  assumed  by  Hasty. 
Note  that  the  previous  model  is  incorporated  into  this  model  as  equations  (3.100), 
(3.101),  (3.102),  (3.105),  and  (3.106).  Yet,  with  the  introduction  of  equations  (3.103) 
and  (3.104),  this  model  now  can  represent  Hasty  et  al’s  mutant  system  describing 
the  lysis  -  lysogeny  decision  of  the  bacteriophage  A  virus. 

The  analysis  of  this  model  is  much  the  same  as  the  previous  model.  It  begins 
with  equations  representing  the  flux  of  one  molecule  through  each  equation. 
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Vi 

=  r 

(3.107) 

V2 

(3.108) 

V3 

=  k.2[DX2]-k2[D\[X2\ 

(3.109) 

V4, 

=  k.3[DX*]-k3[D][X2\ 

(3.110) 

=  k^ilDX^X^]  -  k4[DX2][X2] 

(3.111) 

V6 

=  -krlDX^lP] 

(3.112) 

V7 

=  -kd[X] 

(3.113) 

These  equations  can  then  be  incorporated  into  differential  equations  that  rep¬ 
resent  the  flow  of  each  molecule  through  the  system.  The  simplifled  form  of  this 
system  is  given  below. 


d\X] 

dt 

=  Ui  +  2v2  —  nve  +  U7 

(3.114) 

d[X2] 

dt 

=  ~V2  +  "^3  +  'y4  +  Us 

(3.115) 

dpi 

dt 

=  V3  +  V4 

(3.116) 

d[DX2\ 

dt 

=  -Vs  +  V4 

(3.117) 

d[DX*] 

dt 

=  -Vi 

(3.118) 

d[DX2X2] 

dt 

=  -Vr> 

(3.119) 

substituting  equations 

(3.107)  through  (3.113)  into  equations  (3.114) 

through  (3.119),  the  following  differential  equations  are  created  that  represent  the 
system. 
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®  =  r  +  2(fc_i[X2]  -  h[X]^)  +  nkT[DX2][P]  -  UX] 
^  =  -(fc_i[X2]  -  h[Xf)  +  k.2[DX2]  -  k2[D][X2] 

+  k.s[DX*]  -  /C3p][X2]  +  k.4DX2X2]  -  k4DX2][X2] 
^  =  k_2[DX2]  -  k2[D][X2]  +  k_s[DX*]  -  ks[D][X2] 

=  -{k_2[DX2]  -  k2[D][X2\)  +  k_s[DX;]  -  k4D][X2] 

=  -{k.,[DX*]-k4D][X2]) 

=  -{k.,[DX2X2]-k,[DX2\[X2]) 


With  the  assumption  that  (3.101)  though  (3.104)  represent  fast  reactions, 
steady  state  approximations  can  be  made  regarding  some  of  the  equations  above. 
Specifically,  ^  =  0,  ^  =  0,  =  0,  =  0,  and  =  0.  With  the 

substitutions  x  =  [X],  y  =  [X2],  d  =  [D],  u  =  [DX2],  v  =  [DX2],  z  =  [DX2X2],  and 
Po  =  [P]  the  following  equations  result. 


dx 

dt 

0 

0 

0 

0 

0 


2{k,-iy  —  k^lx)  +  npokru  —  kdX  +  r  (3.120) 

—  (k-iy  —  k^x)  +  k-2U  —  k2dy  +  k-2,v  —  k^dy  +  k-^z  —  k^uy  (3.121) 

k-2U  —  k2dy  +  k-^v  —  k^dy  (3.122) 

—  {k-2U  —  k2dy)  +  k^sv  —  ksdy  (3.123) 

-{k^sv-k^dy)  (3.124) 

—  {k^iZ  —  kiuy)  (3.125) 


Substituting  equation  (3.124)  into  equation  (3.123),  results  in  0  =  —  (/c_2M  ^ 
k2dy).  Then,  substituting  this  result  along  with  equations  (3.123)  and  (3.124)  into 
equation  (3.121)  yields  0  =  —{k^iy  —  kfx)  which  can  be  solved  for  y  as  a  function  of  x 
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so  that  y  =  Then  substituting  y  =  -^x‘^  into  0  =  — (A;_2'U— /c2<^2/),  and  solving 

for  u  yields  the  equation  u  =  Together,  substituting  u  =  '^'^x'^d  and 

y  =  -^x^  into  equation  (3.125)  yields  2:  =  Finally,  substituting 

y  =  ■^x?‘  into  equation  (3.124),  results  in  the  eqiration  v  =  -^-^x'^d.  Thus, 
equations  (3.121)-(3.125)  can  be  written  as  the  following  eqiiations  where  y,  u,  v, 
and  z  are  functions  with  respect  to  the  terms  d  and  x. 


ki  2 

i,  ^ 
k-i 

(3.126) 

ki  k2  2> 
k-i  k—2 

(3.127) 

h  ks  2, 

i,  i,  ^  ^ 

k-i  k-z 

(3.128) 

,  ki  2  ^2  ki  ^ 

(i,  )  i,  i,  ^ 

k—i  k—2  k—4 

(3.129) 

Substituting  these  equations  into  equation  (3.120)  and  simplifying  yields  the 
following  equation. 


dx 

dt 


=  npokr 


h 

k^i 


-^^x^d  —  kdX  +  r 

k^2 


(3.130) 


m-n  =  0,  =  0,  and  =  0  implies  « 


Now,  ™  =  0 


dt 


dlDXj]  ^ 


dt  dt  dt 

d[DX^2^2]  n  4- 1  t r  d[D D D ^ ^2 ^2] 


d[DX2] 

dt 


=  0  or  more  explicitly. 


dt  '  dt  w  VO.  ...VO.V. 

exists  a  constant  number  of  binding  sites  dr  such  that. 


=  0.  This  means  there 


D  +  DX2  +  DX*  +  DX2X2  =  dr 


(3.131) 


which  implies  d -\-  u -\-  v  -\-  z  =  dr 


(3.132) 
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Now,  by  substituting  equations  (3.126)  though  (3.129)  into  equation  (3.132), 
the  following  equation  emerges. 


^1  ^3  2  j  I  /  \2  ^2  ^4  4 


/c_i  /c_3  /c_i  k—2  k—A 


X  d  =  dr 


k_i  /c_2 

Equation  (3.133)  can  then  be  solved  in  terms  of  d  which  results  in 


(3.133) 


d 


dj' 


-x^  + 


/  fci  \2  fc2  kd  4 
k-i  k-3-*'  “'"Vfc-l/  k-2k-i 


(3.134) 


This  equation  can  then  be  incorporated  into  (3.130)  to  give  an  equation  in 
terms  of  x. 


dx 


^1  ^2  ^2  ^4  /^4 

k-2  k-4 


kdx 


(3.135) 


This  last  equation  can  be  simplified  by  the  substitutions  Ki  =  K2  = 


k-2  ’ 


K3  =  7^,  K4  =  7^  and  by  further  substitution  of  K3  =  criA'2,  and  K4  =  a2K2. 


dx 


npQkTKiK2d'rx? 


dt  1  +  (1  +  (Ti)KiK2x‘^  +  a2{Ki)'^{K2)‘^xA 


kdX  +  r 


(3.136) 


Then  using  scaling  techniques  and  the  substitutions  x  =  and  t  =  Tqt, 
this  equation  can  be  rewritten  as 


d^ 


2e2 


npo  krKi  K2dTToXQ  ^ 


kdTo^  + 


Tor 


dr  1  +  (1  +  +  (T2(Ki)2(i^2)2X4C^  '  Xq 


(3.137) 


Then  making  the  substitution  Xq  =  and  Tq  =  the  equation 

further  simplifies. 
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dr 

Finally,  setting  a 
realized. 

_  c^e _ , 

dr  1  +  (1  +  (Ti)^^  +  0-2^"^  ^ 

3.6  Two  New  Models  of  Transcription  and  Translation 

3.6.1  New  Model  #1.  Hasty  et  al.’s  second  model  [13]  is  an  example  of 
how  previous  models  have  been  extended  to  incorporate  different  aspects  of  behav¬ 
ior  in  cellular  systems.  The  model  presented  in  this  section  was  constructed  in  an 
effort  to  represent  transcription  and  translation  in  a  similar  fashion  to  Hasty  et  al.’s 
model [14].  The  goal  of  this  model  is  to  track  the  concentration  of  varying  amounts 
of  two  different  effector  molecules  as  they  relate  to  transcription  and  translation. 
By  incorporating  multiple  binding  of  one  type  of  molecule  into  a  protein  that  effects 
its  own  creation,  an  autoregulatory  feedback  loop  is  created  that  may  produce  inter¬ 
esting  behavior.  In  this  model,  if  m  =  0  and  5  >  0,  then  the  model  is  very  similar 
to  Hasty’s  model  [14].  Although  this  is  a  hypothetical  model,  it  acts  as  a  first  step 
in  a  transition  from  previously  created  models  to  the  construction  of  models  unique 
to  the  Air  Force. 


npQkrpdr 


1  +  (1  +  ryfKfK, 


(3.138) 


_  npokTdr  ^  ^  Hasty  et  al.  eqiiation  is  finally 


3-29 


(3.140) 


aEi  +  Cp^  RNAp 


"  fc_i 

bE2  +  Cr^R  (3.141) 

k-2 

R  +  RNAp  RRNAp  (3.142) 

RNAp  +  DNA  RNAp  +  nEi  +  mi?2  (3.143) 

i?i  degrades  (3.144) 

E2  degrades  (3.145) 


Equations  (3.140),  (3.141),  and  (3.142)  represent  fast  reactions.  In  Eqiiation 
(3.140),  a  is  an  integer  representing  the  number  of  the  particular  type  of  effector 
molecule  Ei  that  bonds  with  the  molecular  complex  Cp  to  form  RNA  polymerase, 
RNAp.  In  Equation  (3.141),  b  is  an  integer  representing  the  number  of  the  particular 
type  of  effector  molecule  E2  that  bonds  with  the  molecular  complex  Cr  to  form 
an  RNA  polymerase  repressor,  R.  Equation  (3.142),  represents  the  bonding  of 
the  repressor  R  to  the  RNA  polymerase  RNAp  to  create  a  complex  RRNAp  that 
prevents  the  RNA  polymerase  from  binding  to  the  DNA  and  activating  transcription 
and  translation. 

Equations  (3.143)  through  (3.145)  are  slow  reactions.  In  Equation  (3.143),  the 
RNA  polymerase  RNAp  attaching  to  the  DNA  promoter  site  to  initiate  transcription 
and  thereby  initiate  translation  which  continues  until  the  RNA  degrades.  The  end 
products  are  the  RNAp,  as  well  as,  the  two  effector  molecules:  Ei,  that  was  created 
n  times  in  translation,  and  E2,  that  was  created  m  times  in  translation.  Equations 
(3.144)  and  (3.145)  represent  the  degradation  of  the  proteins  Ei  and  E2. 

Certain  assumptions  about  the  model  are  made  immediately.  First,  it  is 
assumed  that  the  number  of  molecular  complexes,  Cp  and  Cr,  are  constant.  It  is 
also  assumed  that  the  number  of  DNA  promoter  sites  is  constant.  With  this,  the 
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following  equations  represent  the  flux  of  one  molecule  through  each  of  the  equations 
(3.140)-(3.145). 


=  k-i[RNAp]- hCp[Ei]‘^  (3.146) 

V2  =  k.2[R]  -  k2Cr[E2]'’ 

U3  =  -ks[RNAp][R] 

=  -k^DNA[RNAp] 
vn  =  -^5(^1] 
vq  =  —ke[E2\ 


These  equations  can  be  used  to  formulate  a  system  of  differential  equations 
defining  the  rate  of  change  of  molecular  concentrations  through  the  system.  This 
system  of  differential  equations  is 


d[Ey 

(it 

dm 

dt 

d[RNAp] 

dt 

d[R\ 

dt 


avi  —  nvi  +  U5 


bv2  —  m,Vi  +  vq 


-vi  +vzEVi-Vi  =  -vi  +  U3 


-V2  +  U3. 


(3.147) 


By  substituting  the  system  (3.146)  into  the  system  (3.147),  a  complete  system 
is  generated. 
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dm 

dt 

dm] 

dt 

d{RNAp] 

dt 

d[R] 

dt 


ak_i[RNM  -  akmmT  +  nk^DNA[RNAp]  -  k^mpAAS) 
bk_2[R]  -  bk2Crmf  +  m,kiDNA[RNAp]  -  kem] 

-k_i[RNM  +  hCpmV  -  k3[RNAp]m 

-k^R]  +  ham]''  -  k3{RNAp][R] 


Now,  by  nondimensionalizing,  this  problem  can  be  simplified.  Consider  the 
following  substitutions,  t  =  Tqt,  [Ei\  =  Ei^ei,  [i?2]  =  £-2062,  [RNAp]  =  Pqp,  and 
[R]  =  Ror,  where  |  = 


dm] 

dt 

dm] 

dt 

d[RNAp] 

dt 

d[R] 

dt 


k  TP 

a  ^  V  -  ak^CpToEph^l  +  n 
h/io 

b — - - r  -  bk2CrToE2  €2  +  m, 

£20 


k^DNAToPp 

£io 

kjDNATpPp 

£20 


p  —  /c5Toe(3.149) 
■p  -  kpTpe2 


kiCpToEfr, 

k^.ToP  +  p  -  k3TpRppr 

Ap 

k2CrTpEl^  , 

k-2Tpr  +  -  kpTpPppr 

ito 


To  simplify  system  (3.149),  set  k_2Tp  =  1  so  that  Tp  =  Then  set 

kiCpTpEm  =  1)  which  yields  Ei^  =  With  the  requirement  that 

hCrT„E'^'  =  1,  E,,  =  '-yu  =  Requiring  i=i»  =  1,  yields  P„  = 
Finally,  requiring  —  1,  requires  that  —  r/- 

This  leads  to  a  simplification  of  the  original  system  (3.148). 


3-32 


dei 

(it 

de2 

dt 

dp 

dt 

dr 

dt 


/  ,  hDNA  h 

a[p  -  ej  +  n — - - p  -  - — ei 

k-i  k-2 


b{r  —  63)  +  m 


h£NAC 

pk- 


-p 


k_ 


k- 


1  k\Cpk—\^ 
p-\-  - 


1 

a— 1 


k-2 

hv 


62 


kl. 


2  f''_2 

6-1 


-e?  —  -—pr 

k-2 


-r  + 


k^ 

k-2 


k- 


-pr 


(3.150) 

(3.151) 

(3.152) 

(3.153) 


Now,  because  equations  (3.140),  (3.141),  and  (3.142)  represent  fast  reactions, 
and  RNAp  and  R  only  change  concentrations  within  those  reactions,  they  can  be 
assumed  to  be  at  equilibrium  concentrations.  Therefore,  let  ^  =  0,  and  ^  =  0,  and 
solve  equations  (3.152)  and  (3.153)  for  r  and  p  in  terms  of  ei  and  62  to  get  equations 
that  can  be  plugged  into  equations  (3.150)  and  (3.151).  Therefore,  the  problem  can 
be  formulated  in  terms  of  only  ^  and 


3.6.2  New  Model  ^2.  A  second  problem,  where  only  the  general  derivation 
is  included,  is  listed  here.  This  problem  is  much  the  same  as  the  previous  problem, 
however,  it  separates  transcription  and  translation  into  two  equations.  Transcription 
is  listed  as  equation  (3.154),  and  is  a  slow  equation.  Translation  is  listed  as  equation 

(3.155) ,  and  is  a  relatively  fast  equation  compared  with  the  degradation  of  mRNA 
(denoted  RN Am)  which  is  listed  as  equation  (3. 156)  and  is  a  slow  reaction.  Equation 

(3.155)  is  fast  because  a  large  concentration  of  ribosomes  are  assumed  to  be  in  the 
system. 
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aEi  +  Cp^  RNAp 

k-i 

bE2  +  Cr  ^  R 
k-2 

R  +  RNAp  RRNAp 

RNAp  +  DNA  RNAp  +  RNA^  (3.154) 

RNAjn  +  Ribosomes  RNA^,  +  Ribosomes  +  Ei  +  E2  (3.155) 
RNAm  degrades  (3.156) 

El  degrades 

E2  degrades 

Again,  certain  assumptions  about  the  model  are  made  immediately.  First,  it  is 
assumed  that  the  number  of  molecular  complexes,  Cp  and  CV,  are  constant.  It  is  also 
assumed  that  the  number  of  DNA  promoter  sites  are  constant.  Furthermore,  the 
Ribosomes  are  considered  to  be  a  large  constant.  With  this,  the  following  equations 
represent  the  flux  of  one  molecule  through  each  of  the  previous  equations. 


vi  =  k^ilRNAp]  -  kiCplEi]^  (3.157) 

=  k^2[R]  -  k2Cr[E2]’> 

U3  =  -ks[RNAp][R] 

V4  =  -k4DNA[RNAp] 

U5  =  —k^RibosomeslRN  Am] 

=  -kQ[RNAm] 

V7  =  -krlEi] 
vg  =  -kglE-^ 
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These  equations  can  be  used  to  formulate  a  system  of  differential  equations 
defining  the  rate  of  change  of  molecular  concentrations  through  the  system.  This 
system  of  differential  equations  is 


dm 

dt 

dm 

dt 

d[RNAp] 

dt 


d[R] 

dt 

d[RNAm] 

dt 


avi  -  U5  +  ut 
bV2  -V5  + Vs 

-Ui  +  U3  +  U4  -  U4  =  -Vi  +  V3 
-V2  +  U3 

-U4  +  U5  -  U5  +  ue  =  -U4  +  vq. 


(3.158) 


By  substituting  the  system  (3.157)  into  the  system  (3.158),  a  complete  system 
is  generated. 


dm 

dt 

dm 

dt 

d[RNAp] 

dt 

d[R] 

dt 

d[RNAn^] 

dt 


ak^ddN  Ap]  —  akiCplEiY  +  k3Rihos<ym,es[RN  A^]  —  kj[Ei] 
bkm]  -  bk2Cr[E2t  +  hRibosomes[RNAm]  -  ksm] 

-kmRNAp]  +  kiCp[Ei]^  -  kslRNApWR] 

-k^R]  +  hCm]'"  -  kslRNApWR] 

kiDNA[RNAp\  -  kis[RNA^\ 


Again,  due  to  the  autoregulatory  feedback  of  this  system,  coupled  with  the 
degradation,  this  equation  is  expected  to  show  interesting  behavior.  However,  due 
to  the  complexity  of  this  problem,  it  will  not  be  simplified,  nor  derived  further  here. 
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3. 7  Glucose 


The  glycolytic  pathway  is  metabolic  pathway  that  starts  with  glucose  and 
produces  pyruvate  via  ten  catalyzed  reactions.  “Glycolysis  may  be  considered  to 
occur  in  two  stages.  Stage  I  (Reactions  1-5):  Glucose  is  phosphorylated  and  cleaved 
to  form  two  molecules  of  the  triose  glyceraldehyde-3-phosphate.  This  requires  the 
expenditure  of  two  ATPs  in  an  ‘energy  investment’  (Reactions  1  and  3).  Stage  II 
(Reactions  6-10):  The  two  molecules  of  glyceraldehyde- 3-phosphate  are  converted  to 
pyruvate  with  the  concomitant  generation  of  four  ATPs  (Reaction  7  and  10).”  [34:p 
446]  Figure  3.1  outlines  the  glycolysis  pathway.  Using  both  Gepasi  and  Matlab, 
a  basic  model  of  the  glycolytic  pathway  was  created.  The  model  was  created  by 
using  reaction  equations  associated  with  mass  action  to  define  a  system  of  differential 
equations.  Both  in  Gepasi  and  Matlab,  the  model  is  the  same  so  as  to  compare  the 
Gepasi  simulation  with  the  Matlab  simulation. 

3. 7. 1  Gepasi.  Using  Gepasi  to  create  the  glycolysis  pathway  was  an  easy 
task.  Simply  type  the  twenty  catalytic  reactions  that  define  the  system  into  the 
‘reactions’  box  that  is  listed  on  the  ‘model  definition’  tab.  These  equations  can 
be  written  in  the  form  of  reaction  equations.  Then,  going  to  the  ‘kinetics’  box, 
the  reaction  equations  can  be  described  with  any  one  of  a  series  of  different  kinetic 
types.  The  Glucose  model  used  mass  action.  In  choosing  mass  action,  boxes  appear 
that  allow  for  the  kinetic  constants  to  be  incorporated  into  the  reaction.  Once 
the  reaction  equations  are  complete,  they  are  converted  by  Gepasi  into  differential 
equations  that  can  be  solved  by  the  LSODA  solver  (Livermore  Solver  of  Ordinary 
Differential  Equations).  By  incorporating  initial  conditions  in  the  ‘metabolites’  box 
that  is  listed  on  the  ‘model  definition’  tab,  the  model  is  set  to  run.  By  going  to  the 
‘tasks’  tab,  an  end  time  can  be  chosen  as  well  as  the  number  of  time  steps  LSODA 
will  use  to  approximate  the  system.  By  checking  the  first  two  boxes  in  the  report 
section  of  the  ‘tasks’  tab,  the  output  data  from  Gepasi  will  include  the  concentrations 
of  all  reactant  species  at  the  end  time  chosen. 
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ATP  ADP 


Mg++,  HexQkinase 


Beta-D-Fructose-1,  fi-Bisphosphate 


ATP 

i  3 


ADP 

1 


2  I  Phosphoglucoisomerase 


Mg++,  Phosphofhictokinase 
Aldolase 


B  eta-D  -Fructos  e-  d  -Pho  sphate 


Gly  c  erald  ehy  de-  3  -Phosphate 


Glyc  erald  ehyde-  3  -Pho  sphate 
Dehydrogenase  g 


1 , 3  -Bispho  spho  glyc  erate 


Tfiose  Phosphate  Isomerasi 
2  (NAD-h,  Pi) 

2  (nADh) 

Mg++,  Pho  spho  glyc  erate  Kinase 


Dihydroxyacetone  Phosphate 


t  ^  i 

2  ADP  2  ATP 


3 -Pho  spho  glyc  erate 


g  Pho  spho  glyc  erate  Mutase 


Mg-r-r,  P3mjvate  Kinase 


Mg-l-^,  Enolase 


Figure  3.1.  The  glycolysis  metabolic  pathway 
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3. 7.2  MATLAB.  In  the  Matlab  script,  there  are  36  different  molecules  that 
are  changing  in  the  system.  Each  reaction  has  a  reaction  rate  coefficient  associated 
with  it.  The  reaction  rate  coefficients  are  given  in  the  code  by  kl  —  k37.  The  forward 
reaction  rate  coefficients  are  3  and  the  backwards  reaction  rate  coefficients  are  1.  In 
this  model,  all  reaction  rate  coefficients  are  fixed.  In  futiire  models,  they  may  be 
considered  variables  or  functions.  Part  of  the  script  is  listed  below.  In  reaction  ^2 
there  is  only  one  kinetic  constant  because  there  is  no  backward  reaction,  the  reaction 
is  one-way. 


kl  =  3; 

REACTION#! 

(3.159) 

k2  =  1; 

(3.160) 

k3  =  3; 

REACTION#2 

(3.161) 

kA  =  3; 

REACTION#3 

(3.162) 

=  1; 

(3.163) 

Each  reaction  is  dependent  on  the  molecules  formed  by  previous  reactions. 
Therefore,  all  reaction  in  this  system  are  interconnected  through  time.  This  inter¬ 
connection  makes  modeling  the  system  as  a  whole,  a  very  difficult  process.  How¬ 
ever,  first  looking  at  each  reaction  individually,  then  incorporating  the  equations 
for  the  individual  reactions  into  the  complete  model  makes  this  model  easier  to  con¬ 
struct.  There  are  10  individual  reactions  in  the  very  simple  biochemical  model  of  the 
metabolic  pathway,  Glycolysis.  However,  each  reaction  is  a  catalytic  reaction.  This 
being  the  case,  each  reaction  can  be  split  into  two  reactions  that  more  accurately 
describe  the  catalytic  reaction.  [18]  This  changes  the  10  reactions  into  20  slightly 
more  complex  reactions  involving  hypothetical  molecular  states  of  catalytic  transi¬ 
tion.  Each  of  the  20  reactions  has  a  reaction  rate  depending  on  the  concentration 
of  molecules  entering  into  the  system  and  the  concentration  of  those  leaving.  Before 
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looking  at  the  system  as  a  whole,  the  individual  reaction  eqiiations  are  constructed. 
Below  is  a  summary  of  all  the  reactants  of  the  system.  Also  included  are  the  first 
two  of  twenty  reaction  equation  along  with  the  biochemical  eqiiations  that  explain 
them. 

_ QUICK  KEY _ 

ENZYMES  =  1/(1)  -  1/(9),  //(35) 

CATALYST  TRANSITION  =  y{10)  -  y{19) 

ENERGY  =  //(20)  -  y(24) 

SUBSTRATES  =  y(25)  -  y(34),  y{36) 


The  first  two  reaction  equations  in  the  system  explain  the  reaction  process  that 
takes  alpha-D- Glucose  to  alpha-D-Glucose-6-phosphate  via  the  catalyst  Hexokinase 
and  coenzyme  Mg+^.  In  the  process,  this  reaction  requires  energy  in  the  form  of 
ATP  and  creates  the  by-product  ADR 

The  first  reaction  equation  takes  alpha-D-Glucose  in  the  presence  of  ATP, 
Hexokinase,  and  Mg++  and  produces  Gatalyst_transition_l.  It  is  presented  as 
follows. 

alpha-D-Glucose  +  ATP  +  Mg^^+  Hexokinase  {k2)  <==>  {kl)  Cata- 
lyst_transition_l 

The  algebraic  equation  corresponding  to  it  is 

egl  =  —kl  *  y{25)  *  y{20)  *  y{l)  *  y{2)  +  k2  *  y(10). 

The  second  reaction  equation  describes  Catalyst_transition_l  producing  alpha- 
D-Glucose-6-phosphate,  ADP,  Hexokinase,  and  Mg++.  It  is  presented  as  follows. 
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Catalyst_transition_l  =>  (^3)  alpha-D-Glucose-6- phosphate  +  ADP  + 
Hexokinase 

The  algebraic  equation  corresponding  to  it  is 

eg2  =  — A;3  *  2/(10). 

The  concentration  of  each  of  the  36  molecules  change  with  respect  to  time. 
Some  of  these  molecules  exist  in  more  than  one  equation.  Each  of  the  following 
equations  is  a  reaction  rate  written  with  respect  to  a  single  moleciile  on  the  left 
hand  side  of  the  biochemical  equation.  The  negative  of  that  eqiiation  will  give  the 
reaction  rate  written  with  respect  to  a  single  molecule  on  the  right  hand  side  of  the 
biochemical  equation.  To  find  the  reaction  rate  of  a  specific  molecule  over  the  entire 
system,  add  each  equation  where  the  molecule  appears  on  the  left  hand  side  of  the 
biochemical  equation.  Of  course,  if  the  molecule  appears  on  the  right  hand  side  of 
the  biochemical  equation,  add  the  negative  of  the  above  equations. 


dyjl) 

(it 

dy{2) 

(it 


eql  —  eq2  +  eq5  —  eqQ  +  eql3  —  eqlA  +  eqll  —  egl8  +  egl9  —  eg20; 
egl  -  eg2;  (3.164) 


The  first  equation  (3.164)  from  the  larger  system  exhibits  a  change  of  Mg++ 
with  respect  to  time.  The  second  equation  (3.164)  from  exhibits  a  change  of  Hex¬ 
okinase  with  respect  to  time.  Looking  through  this  system,  you  will  notice  that 
these  equations  are  nonlinear.  This  makes  solving  the  system  exactly,  virtually  im¬ 
possible.  However,  a  numerical  solver  such  as  Matlab’s  ‘ode23’  or  ‘ode23s’  can  solve 
the  system  easily.  The  script  written  for  Glucose  uses  the  ‘ode23s’  stiff  solver  to 
solve  the  system.  This  stiff  solver  requires  a  set  of  initial  conditions  for  all  thirty  six 
reactant  species  as  well  as  a  time  interval  in  order  to  solve  the  system.  Comparing 
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the  output  of  the  Matlab  script  at  a  specific  end  time  to  the  output  of  Gepasi,  the 
two  systems  respond  the  same.  Outpnt  from  Gepasi  and  Matlab  is  listed  in  chapter 

4. 

3. 8  Summary 

This  chapter  looked  at  previonsly  developed  methods  for  determining  the 
steady  state  behavior  for  systems  of  antonomons  ordinary  differential  eqiiations. 
These  methods  included  defining  stability  reqnirements  for  linear  models  as  well  as 
linearizing  nonlinear  models  to  test  for  stability  in  local  regions.  It  then  used  these 
methods  to  analyze  chemical  systems,  snch  as  the  Brusselator  and  Schnackenberg’s 
model,  by  reformulating  the  chemical  models  into  systems  of  differential  eqiiations. 
A  stability  analysis  was  performed  on  these  stable  states  to  determine  the  sensitivity 
of  the  steady  state  behavior  under  different  parameter  settings.  By  considering  in¬ 
teresting  chemical  systems  with  bifurcations,  a  variety  of  different  types  of  behavior 
were  described.  Finally,  by  evaluating  models  of  transcription  and  translation,  it  has 
been  shown  that  these  interesting  reaction  models  can  be  incorporated  into  biological 
systems  to  define  specific  aspects  of  biological  behavior.  With  a  better  understand¬ 
ing  of  different  forms  of  stability,  such  as  asymptotic  and  oscillatory  stability,  and 
different  techniques  for  finding  that  stability,  to  include  linearization  techniques  and 
substitutions,  it  will  be  easier  to  construct  models  that  when  incorporated  together, 
represent  real  biological  systems  that  can  be  tested  to  fulfill  Air  Force  needs. 
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IV.  Model  Results 


4-1  OvervieAV 

A  number  of  the  models  constructed  in  Chapter  3,  the  analysis  portion  of 
this  thesis,  potentially  exhibit  interesting  behavior  siich  as  limit  cycle  solutions  and 
multiple  stable  steady  state  solutions.  Under  specific  conditions,  this  behavior  may 
be  sensitive  to  variable  concentrations  or  model  parameters.  This  Chapter  will 
investigate  changes  of  the  small  linear  system  model  of  Chapter  3,  the  Brusselator 
model  [27],  the  Schnackenberg  model  [30],  the  model  constriicted  by  Hasty  et  al. 
[14],  a  reaction  inside  the  Glucose  metabolic  pathway,  and  the  larger  Glucose  model 
[34]  as  a  result  of  perturbations  of  initial  conditions  and  parameters.  In  this  chapter, 
all  quantities  are  nondimensional. 

4-2  A  small  linear  system 

The  small  system  constructed  in  Chapter  3,  namely 

^  =  -(3  +  £)A  +  3F  (4.1) 

dY  _  (2  +  £)\  ^ 

was  analyzed  to  find  the  values  of  the  perturbation  term  e  where  steady  state  be¬ 
havioral  changes  take  place.  In  this  Chapter,  a  numerical  analysis  of  specific  cases 
adds  more  detail  to  how  the  concentrations  of  X  and  Y  behave  given  different  values 
for  e.  The  basic  analysis  of  this  linear  system  leads  to  a  greater  understanding  of 
how  a  linearized  system  will  behave  in  a  neighborhood  of  a  steady  state. 

The  analysis  of  Chapter  3  shows  that  any  value  e  belonging  to  the  interval 
(^8/3,0)  will  lead  to  a  stable  solution.  Figure  4.1  represents  X  and  Y  when  e  is 
given  the  value  e  =  —1,  and  the  initial  conditions  for  X  and  Y  are  Xq  =  To  =  5. 
The  eigenvalues  of  the  system  under  these  conditions  are  Ai  =  —.5,  and  A2  =  ^2.5 
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Figure  4.1.  Steady  state  stability  in  a  linear  system  when  £  =  —  1 

which  implies  the  solution  is  stable  as  predicted  in  Chapter  3.  With  50  time  steps, 
the  final  concentration  values  for  X  and  Y  reach  the  steady  state  value  of  X  =  0 
and  y  =  0  asymptotically. 

The  analysis  of  Chapter  3  also  shows  that  any  positive  value  for  e  will  lead 
to  an  unstable  solution.  Figure  4.2  represents  X  and  Y  when  e  is  given  the  value 
£  =  1,  and  the  initial  conditions  for  X  and  Y  are  Xq  =  =  5.  The  eigenvalues  of 

the  system  are  Ai  =  .5,  and  A2  =  —5.5  which  implies  the  solution  is  a  saddle.  With 
5  time  steps,  the  final  concentration  values  for  X  and  Y  moves  quickly  away  from 
the  steady  state  value  to  a  concentration  of  X  =  45.8812  and  Y  =  68.8218. 

This  saddle  point  behavior  is  seen  again  in  when  £<-§.  To  illustrate  this,  £ 
was  given  the  value  £  =  —2.9,  with  the  initial  conditions  for  X  and  Y  left  at  Xq  = 
Iq  =  5.  The  eigenvalues  under  these  conditions  are  Ai  =  —1.45,  and  A2  =  0.35 
which  again  implies  the  solution  is  a  saddle  as  predicted  in  Chapter  3.  With  5  time 
steps,  the  final  concentration  values  for  X  and  Y  moves  away  from  the  steady  state 


4-2 


70 


Figure  4.2.  Steady  state  instability  in  a  linear  system  when  e  =  1 

value  to  a  concentration  of  X  =  69.6323  and  Y  =  10.4475.  This  is  plotted  in  Figure 
4.3. 

These  models  show  both  a  case  of  asymptotic  stability  to  a  steady  state,  and 
two  cases  of  divergent  behavior  in  the  form  of  a  saddle  point.  They  also  indicate  that 
solutions  of  a  system  can  be  highly  sensitive  to  small  perturbations  of  parameters. 
This  analysis  gives  a  partial  description  of  the  behavioral  characteristics  of  a  linear 
system.  This  analysis  can  then  be  expanded  to  evaluate  more  complex  systems 
through  the  process  of  linearization. 

4-3  The  Brusselator 

In  equations  (4.2),  the  differential  equations  representing  the  Brusselator  equa¬ 
tion  include  parameters  that  can  be  perturbed  just  as  the  previous  linear  model. 
These  parameters  in  the  Brusselator,  represented  as  a  and  h,  are  controlled  values 
that,  when  perturbed,  can  introduce  very  interesting  behavior  to  the  system.  This 
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Figure  4.3.  Steady  state  instability  in  a  linear  system  when  e  =  —2.9 

section  will  analyze  the  behavior  of  the  Brusselator  to  distinguish  the  difference  in 
behavior  of  a  nonlinear  systems  from  that  of  linear  systems. 


dx 

dt 

dt 


a  —  (6  +  l)x  +  x'^y 
bx  —  x?y 


(4.2) 


In  Figure  4.4,  a  =  0.60710,  and  b  =  1.5,  so  the  eigenvalues,  based  on  the 
linearization  of  section  3.3.2,  are  Ai^2  =  6.5711  x  10“^  ±  .603  54^  which  implies 
cyclic  divergence  from  the  steady  state  point  located  at  {x,y)  =  (0.6071,2.470  7). 
However,  once  outside  a  small  neighborhood  of  the  steady  state,  the  concentrations 
approach  a  limit  cycle.  This  is  not  unexpected  since  linearization  indicates  only 
local  behavior. 

As  another  example,  consider  the  case,  a  =  0.8071  and  b  =  1.5,  so  the  eigen¬ 
values  are  Ai,2  =  —7.  5705  x  10^^  ±  .  803  54i,  which  implies  the  steady  state  located 
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Figure  4.4.  Limit  cycle  behavior  in  the  Brusselator 

at  (x,  y)  =  (0.8071, 1.8585)  is  an  oscillatory  stability  point.  In  Figure  4.5,  the  initial 
concentrations  are  located  at  (xo,yo)  =  (0.9071, 1.9585),  and  are  therefore  within  a 
sufficiently  small  neighborhood  of  the  steady  state  point  to  allow  linearization  to  be 
used  as  a  means  of  determining  the  concentrations  behavior.  As  the  concentrations 
change  with  time,  they  approach  the  steady  state. 

Although  the  Brusselator  model  is  a  hypothetical  chemical  oscillator,  the  be¬ 
havior  within  it  may  exist  inside  larger  systems  that  are  models  of  cellular  functions. 
By  understanding  how  linearity  affects  this  system,  more  awareness  can  be  gained 
about  the  conditions  linearity  imposes  on  larger  systems.  This  may  eventually  lead 
to  a  better  understanding  of  large  models  involving  cellular  processes. 

4-4  Schnackenberg 

The  Schnackenberg  model,  presented  in  Chapters  2  and  3,  has  stability  charac¬ 
teristics  that  are  important  to  understand  when  developing  large  biological  models. 
This  section  will  do  a  sensitivity  analysis  on  both  a  limit  cycle  and  stable  steady 
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Figure  4.5.  Cyclic  stability  in  the  Brusselator 

state  within  the  Schnackenberg  model.  The  sensitivity  analysis  will  consist  of  per¬ 
turbations  of  reactant  concentrations  from  stable  steady  state  and  limit  cycle  as  well 
as  perturbations  of  constant  parameters  within  the  model.  Each  sensitivity  analysis 
is  evaluated  by  a  plot  showing  how  concentrations  change  through  time.  In  these 
plots,  the  concentration  of  X  is  measured  on  the  x-axis,  and  the  concentration  of  Y 
is  measured  on  the  y-axis.  This  analysis  will  lead  to  understand  how  systems  may 
be  very  sensitive  to  small  perturbations. 

Figure  4.6  indicates  limit  cycle  behavior  for  the  Schnackenberg  model  given  by 
equation  (3.49)  and  (3.50).  This  model’s  constant  parameter  values  are  a  =  0.2091 
and  b  =  .1.  From  equations  (3.53)  and  (3.54)  the  model  has  a  steady  state  located 
at  the  point  {X,Y)  =  (0.3091,2.1884).  From  the  characteristic  equation  (3.56), 
the  eigenvalues  have  a  positive  real  part,  so  the  steady  state  is  unstable.  The 
initial  condition  for  Figure  4.6  is  (Yq,  Yq)  =  (0.4091,  2.2884).  This  leaves  the  initial 
condition  within  the  interior  of  the  limit  cycle  which  is  indicated  in  Figure  4.6.  If 
the  initial  point  is  any  point  inside  the  limit  cycle,  other  than  the  steady  state  point, 
the  resulting  trajectory  will  trace  a  path  out  to  the  limit  cycle. 
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Figure  4.6.  Interior  convergence  to  a  limit  cycle  in  a  Schnackenberg  model 

Figure  4.7  shows  the  same  limit  cycle  with  initial  condition  (Xq,  Fq)  =  (1-3091,  3.1884). 
This  puts  the  initial  condition  outside  the  limit  cycle.  Note  that  the  trajectory  which 
the  solution  traces  moves  away  from  the  limit  cycle  before  it  moves  to  it. 

Figure  4.8  shows  another  type  of  behavior  for  the  Schnackenberg  model. 

This  model’s  constant  parameter  values  a  =  0.1091  and  b  =  0.1.  From  equa¬ 
tions  (3.53)  and  (3.54)  the  model  has  a  steady  state  located  at  the  point  (X,  T) 

=  (0.2091,2.4952).  From  equation  (3.56),  the  eigenvalues  have  negative  real  parts 
so  the  steady  state  is  stable.  The  initial  condition  in  this  model  is  (Xo,yo)  = 
(0.3091,2.5952).  This  leaves  the  initial  condition  only  a  small  distance  from  the 
stable  steady  state  point.  The  concentration  moves  much  further  from  the  steady 
state  before  returning  back.  As  the  concentration  moves  back  towards  steady  state, 
its  trajectory  wraps  tightly  around  the  steady  state  point.  To  graphically  observe 
convergence  to  the  steady  state  point,  every  tenth  and  then  every  50th  point  was 
plotted. 
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Figure  4.7.  Exterior  convergence  to  a  limit  cycle  in  a  Schnackenberg  model 

In  Figure  4.8,  the  graph  shows  that  the  system  is  very  sensitive  to  small  per¬ 
turbations.  In  a  biological  model,  this  type  of  behavior  might  be  detrimental  to  a 
cell.  In  Figure  4.8,  the  small  perturbation  leads  the  concentrations  of  the  reactants 
to  toxic  levels  that  may  kill  or  damage  the  organism.  In  this  case,  it  becomes  very 
important  for  the  larger  biological  model  to  regulate  this  process  so  as  to  prevent  its 
functions  from  loosing  control.  On  the  other  hand,  this  type  of  behavior  may  be 
integrated  into  a  cell  model  as  an  amplifier  that  triggers  a  switch  in  later  reactants 
from  one  steady  state  concentration  to  another.  This  type  of  behavior  is  described 
in  the  next  section. 

The  final  Schnackenberg  example  shows  stable  behavior  (Figure  4.9).  This 
model’s  constant  parameter  values  are  a  =  0.0091  and  h  =  0.1.  The  model  has 
a  stable  steady  state  located  at  the  point  {X,Y)  =  (0.1091,0.7679).  The  initial 
condition  in  this  model  is  (Xo,To)  =  (0.2091,0.8679). 

The  stability  in  this  example  is  asymptotic  which  is  quite  different  from  the 
oscillatory  convergence  or  the  convergence  to  a  limit  cycle  of  the  previous  examples. 
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Figure  4.8.  Oscillator  convergence  in  the  Schnackenberg  model 

In  this  and  the  previous  examples,  the  steady  state  point  is  based  on  the  constant 
parameter  values  using  equations  (3.53)  and  (3.54).  Furthermore,  each  example  had 
an  initial  condition  placed  at  (Xo,Fo)  =  (-^ss  +  SiiF'ss  +  £2)  where  ei  =  £2  =  0.1. 
The  only  perturbation  in  each  model  came  from  the  model’s  constant  parameter 
a.  Therefore,  the  behavior  of  the  Schnackenberg  model  is  obviously  sensitive  to 
perturbations  from  the  constant  parameter  values. 

The  large  behavioral  changes  due  to  small  perturbations  of  parameters  in  the 
small  Schnackenberg  model  give  insight  into  the  great  complexity  of  larger  systems. 
From  the  Schnackenberg  model,  it  becomes  very  clear  that  larger  biological  systems 
need  to  control  the  range  of  perturbation  on  their  parameter  values  so  as  to  prevent 
extreme  changes  in  the  system  that  could  kill  the  organism. 

4.5  Hasty  et  al.  ’s  Small  Model 

Hasty  et  al.’s  model,  developed  in  Chapters  two  and  three,  was  an  example  of 
a  bistable  system.  In  this  section.  Figures  4.10  and  4.11,  representing  a  stability 
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Figure  4.9.  Asmptotic  convergence  in  the  Schnackenberg  model 

analysis,  will  give  insight  into  exactly  how  the  system  behaves  under  small  perturba¬ 
tions  in  the  concentrations  of  reactants  X  and  X2  from  their  steady  state.  Figures 
4.12  and  4.13,  representing  bifurcation  diagrams,  will  show  exactly  how  the  steady 
state  concentration  of  the  reactant  X  behaves  under  different  values  of  the  constant 
parameters  a  and  S.  Using  both  sensitivity  graphs  and  bifurcation  diagrams,  the 
behavior  of  Hasty  et  al.’s  model  will  be  fully  understood. 

Figure  4.10  shows  stable  behavior.  Its  constant  parameter  values  are  a  =  50, 
(5  =  0.4,  Ki  =  1,  K2  =  1,  13  =  K1K2,  7  =  10,  and  dT  =  10.  For  these  parame¬ 
ter  values,  the  model  has  a  three  steady  states  located  at  the  points  {X,X2)i  = 
(0.0552,0.0030),  {X,X2)2  =  (0.1499,0.0225),  (A, ^2)3  =  (4.8349,23.3761)  where  X 
is  respectively  one  of  the  roots  of  equation  (3.99)  and  X2  satisfies  equation  (3.91)  with 
^  =  Ki  =  1.  The  initial  condition  in  this  model  is  (X,  ^2)0  =  (0.0652,0.0130). 
Consequently,  the  initial  condition  is  a  small  distance  from  the  steady  state  point 
(X,  X2)i.  As  is  observed  in  Figure  4.10,  the  concentrations  quickly  return  to  the 
original  steady  state.  Thus,  it  is  clear  that  the  steady  state  at  (X,  X2)i  is  stable. 
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Figure  4.10.  Lower  steady  state  convergence  in  the  Hasty  model 

Figure  4.11  also  shows  stable  behavior.  All  parameter  values  are  the  same 
as  for  the  previous  example.  The  only  change  is  that  (A,  ^2)0  =  (0.1552,0.1030). 
In  fact,  this  puts  steady  state  (A,  ^2)2  between  (A,  A2)i  and  initial  condition.  As 
is  observed  in  Figure  4.11,  the  concentration  moves  away  from  both  steady  state 
(A,  A2)i  and  steady  state  (A,  A2)2,  and  converges  to  steady  state  (A,  A2)3.  This 
indicates  that  the  steady  state  at  (A,  A2)2  is  unstable,  and  steady  state  (A,  A2)3  is 
stable.  Thus,  these  numerical  examples  are  consistent  with  the  analysis  by  Hasty 
et  al.  concluding  that  this  is  a  bistable  system. 

The  next  two  hgures  are  bifurcation  diagrams  that  explain  the  behavior  of  the 
steady  state  concentrations  as  a  function  of  a  single  parameter.  Figure  4.12  describes 
the  steady  state  concentration  of  X  as  a  function  of  a.  The  constant  parameter 
values  are  8  =  0.4,  Ai  =  1,  A2  =  1,  /I  =  K1K2,  7  =  10,  and  dT  =  10.  For  each  a,  the 
values  of  [A]  are  determined  by  solving  equation  (3.99)  where  u  =  [A].  The  constant 
parameter  a  may  change  for  a  number  of  reasons.  Recall  from  equations  (3.95)  and 
(3.96)  that  a  =  nksPoPdT,  so  if  the  concentration  of  RNA  polymerase  po  were  to 
somehow  change,  a  would  change  in  proportion  to  it.  As  indicated  previously,  when 
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Figure  4.11.  Upper  steady  state  convergence  in  the  Hasty  model 

a  =  50,  the  steady  state  points  are  at  approximately  {X,X2)i  =  (0.0552,0.0030), 
(X,  ^2)2  =  (0.1499, 0.0225),  and  {X,  ^2)3  =  (4.8349, 23.3761).  Also,  by  the  stability 
analysis,  it  has  been  show  which  steady  states  are  stable  and  which  are  unstable. 
The  graph  shows  that  bifurcation  points  exist  at  approximately  a  =  20,  and  a  =  63. 
At  a  =  20,  two  of  the  steady  state  solutions  go  from  complex  to  real  values  and  then 
at  a  =  63  two  steady  states  go  from  real  to  complex  values.  Thus  when  a  <  20, 
there  is  one  real  steady  state  value,  when  20  <  a  <  63,  there  are  three  real  steady 
state  values,  and  again  when  63  <  a,  there  is  one  real  steady  state  value.  Figure 
4.12  is  exactly  the  same  as  Figure  la  of  Hasty  et  al. 

This  bifurcation  plot  can  also  be  associated  with  8.  Recall  from  equation 
(3.74)  that  8  is  the  basal  rate  of  production  of  X.  This  may  change  depending  on 
occurrences  outside  of  the  system.  For  8  varying  about  the  x-axis,  a  =  25,  and  all 
other  parameters  held  constant  with  rates  described  previously,  it  is  obvious  that 
small  changes  in  8  lead  to  large  changes  in  the  system.  For  example,  at  (5  =  .4,  there 
are  three  stable  states.  Yet,  by  changing  8  such  that  (5  =  1.4,  there  is  only  one  real 
steady  state.  This  implies  that  5  is  a  very  sensitive  parameter. 
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Figure  4.12. 


Figure  4.13. 


A  Hasty  model  bifurcation  plot  that  maps  the  steady  state  of  [X]  with 
respect  to  a 


A  Hasty  model  bifurcation  plot  that  maps  the  steady  state  of  [X]  with 
respect  to  6 
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Bistability  in  Hasty  et  al.’s  model  shows  how  perturbations  can  permanently 
affect  a  system.  Furthermore,  through  stability  analysis,  the  behavior  of  each  steady 
state  can  be  described,  and  through  bifurcation  diagrams,  these  steady  states  can 
be  tracked  with  different  concentrations  of  parameters.  In  larger  models,  bistability 
may  act  as  a  regulatory  switch,  to  control  the  model’s  behavior. 

4.6  Glucose 

4.6.1  Small  Glucose  Model.  An  evaluation  on  a  small  subsystem  within 
the  Glucose  model  showed  interesting  behavior.  The  subsystem  evaluated  was  the 
reaction  taking  a-D-Glucose-6-phosphate  to  /?-D-Fructose-6-phospate  in  the  pres¬ 
ence  of  Phosphoglucoisomerase  or  the  reaction  taking  3-Phosphoglycerate  to  2- 
Phosphoglycerate  in  the  presence  of  Phosphoglycerate  Mutase.  The  reaction  equa¬ 
tions  for  the  first  small  reaction  system  are 


a-D-Glucose-6-phosphate  (yi)  +  Phosphoglucoisomerase  (^2)  ^ 

kf, 

Catalyst  Transition  2  (ya) 


Catalyst  Transition  2  (ys) 

k7 

/3-D-Fructose-6-phospate  {y^)  +  Phosphoglucoisomerase  (y2). 


This  leads  to  differential  equations  of  the  form 


dyi 

dt 

dy2 

dt 

dm 

dt 

dm 

dt 


-hym  +  hys 

-kmy2  +  (^5  +  ^6)2/3  -  k7y2y4 
kmy2  -  (^5  +  ^6)2/3  +  ^72/22/4 
keys  -  ^72/22/4 


(4.3) 
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since,  ^  +  ^  =  0,  ^2  +  2/3  =  c  where  c  =  ^2(0)  +2/3(0).  Thus,  2/3  =  0^2/2- 

Substituting  into  (4.3)  results  in  the  equations 


dyi 

dt 

(kn 

dt 

dy4 

dt 


-k^yiy2  +  h{.c  -  2/2) 

-^42/12/2  +  [h  +  h){c  -  2/2)  -  /i;72/22/4 
^6(c-  2/2)  -  ^7 2/22/4- 


(4.4) 

(4.5) 

(4.6) 


Steady  state  values  can  be  located  by  setting  ^  =  0,  ^  =  0,  and  ^  =  0. 
The  resulting  equations  are 


0  = 

-hym  +  h{c  -  y2) 

(4.7) 

0  = 

-^42/12/2  +  {h  +  h){c  -  2/2)  -  k7y2y4 

(4.8) 

0  = 

k&{c  -  2/2)  -  ^72/22/4 

(4.9) 

Let  2/2  =  a  and  substitute  into  equations  (4.7)-(4.9). 


0  = 

—k^yia  +  kr,{c  —  a) 

(4.10) 

0  = 

—k^yia  +  {k^  +  kQ){c  —  a)  —  k^ay^ 

(4.11) 

0  = 

k(i{c  —  a)  —  k^ay^ 

(4.12) 

Solve  equations  (4.10)  and  (4.12)  for  2/1  and  2/4  to  get  steady  state  concentrations  in 
terms  of  a. 
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(4.13) 


yi 

2/4 


h{c-  a) 
k^a 

k(s{c- a) 
k^a 


(4.14) 


Equation  (4.11)  is  a  linear  combination  of  equations  (4.10)  and  (4.12)  so  a  remains 
an  arbitrary  value.  Therefore,  the  steady  state  values  for  yi,y2  ,2/3  ,  and  2/4  are  given 
by 


2/1 

kf,{c  -  a) 
k^a 

(4,15) 

2/2 

=  a 

(4.16) 

2/3 

=  c  —  a 

(4.17) 

2/4 

kQ(c  —  a) 
kja 

(4.18) 

Figure  4.14  describes  the  steady  state  concentration  of  a-D- Glucose-6- phosphate 
(2/1),  and  /3-D-Fructose-6-phospate  (2/4)  as  a  function  of  the  steady  state  concentra¬ 
tion  of  Phosphoglucoisomerase  (2/2)-  In  this  graph,  c  =  100,  k^  =  2,  k^  =  1,  ke  =  2, 
kj  =  1,  and  a  lies  within  the  range  [25,  75].  As  values  for  a  change,  a  steady 
state  line  is  created  that  defines  the  steady  state  concentrations  for  the  system.  As 
observed  from  the  equations  (4. 15)- (4. 18),  different  values  for  the  kinetic  constants 
lead  to  different  shapes  for  the  steady  state  curve.  In  the  next  section,  the  small 
model  will  be  incorporated  into  a  large  model  of  the  glucose  pathway.  An  example 
which  is  consistent  with  incorporating  the  small  model  into  the  large  model  includes 
initial  conditions  2/1  =  50,  2/2  =  100,  2/3  =  0,  and  2/4  =  0.  This  yields  a  steady  state 
value  of  2/1  =  0.456,  2/2  =  52.281,  2/3  =  47.  718,  and  2/4  =  1.825. 

Testing  the  steady  state  behavior  of  the  steady  state  curve  under  different 
values  of  the  kinetic  constants  showed  that  the  steady  state  curve  was  stable.  That 
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Figure  4.14.  Steady  state  curve  for  the  small  glucose  model 

is,  values  that  were  pushed  from  the  steady  state  line,  returned  to  it.  Unfortunately, 
tracking  the  system  to  a  steady  state  point  is  more  difficult  than  simply  substituting 
in  an  initial  concentration  a  for  Phosphoglucoisomerase  (^2)-  If  a  concentration 
within  the  system  is  moved  away  from  its  steady  state  point  on  the  steady  state 
curve,  then,  in  the  process  of  moving  back  to  the  steady  state  curve,  a  changes. 
Thus,  the  system  will  move  to  a  new  point  on  the  steady  state  curve.  So,  the 
system  is  stable  with  respect  to  the  steady  state  curve,  but  not  with  respect  to 
individual  steady  state  points.  This  is  indicated  in  Figure  4.15. 

Further  analysis  on  the  small  system  determined  whether  there  was  different 
types  of  behavior  near  the  steady  state  for  different  combinations  of  positive  real 
kinetic  constants. 

First,  the  system  (4.3)  was  simplified  to  (4.4)-(4.6)  where  the  third  equation  of 
system  (4.3)  was  omitted  because  j/3  =  c—yi-  The  Jacobian  of  the  system  (4.4)-(4.6) 
takes  the  form 
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Figure  4.15.  Convergence  to  a  new  point  on  the  steady  state  curve  in  the  small 
glucose  model 


^  -hy2  -hVi  -h  0  ^ 

•/  =  -hy2  -hyi  -  (h  +  h)  -  hVi  -hVi 

y  0  -^6  -  ^72/4  -kyXj4^ 

Then  solving  for  the  characteristic  equation,  det|A/  —  J|  =  0 


det 


^  A  +  k4tj2  k^yi  Th  0  ^ 

km  A  +  km  +  {h  +  ko)  +  km  km 
0  kQ-\-  km  A  +  kry^  j 


=  0 


yields  the  cubic  equation 


(4.19) 


(4.20) 


A^  +  (2k’jy4  +  k^yi  +  ^5  +  +  k4y2)X^  +  OqA  —  0. 

where  Oq  =  k^yik^y^  +  k^km  +  “^kmkm  +  kmk^  which  simplifies  to  a  quadratic 
by  factoring  out  a  A.  Thus  the  first  root  is 
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A  =  0 


(4.21) 


and  the  quadratic  equation  is 


A  +  {‘2k, +  k^yi  k^  kg  A;4|/2)A  +  oo  —  0  (4.22) 

where  oq  =  /c42/i/c72/4  +  /c5/i^72/4  +  2kAy2hyA  +  kAy2kQ. 

The  solution  to  the  quadratic  takes  the  form  A  =  —k^y^  —  ^k^yi  —  \k^  —  \kQ  — 
\kAy2±\'/Di,  where  Di  =  2klyiy2+4:k7yAke-4:k4y2k7yA+4:k^yl-2kAy2k6+2kAyik6  + 
2/c4yi/c5  +  2/c5/c4y2  +  kjyf  +  kf  +  2k^k(i  +  kl  +  kjyl 

Substituting  the  steady  state  points  in  for  yi,  y2,  and  7/4  given  by  eqiiations 
(4.15),  (4.16),  and  (4.18)  and  simplifying  gives  the  value  A  =  — ■^(2A:6C—  kQa  +  k5C  + 
k^a'^)  ±  where  ri  =  2k4k^a‘^c  +  kla^  +  2k4a^kQ  +  d/sgC^  +  2k^kQac  +  + 

/cla^  —  4/cgQ!C  —  Ak^k^a^c 

If  ri  <  0,  then  the  system  has  oscillatory  behavior.  Now,  ri  can  be  simplified 
by  omitting  some  positive  terms.  This  can  be  accomplished  by  setting  k^  =  0.  The 
new  equation  is 


/cga^  +  2kikQa^  +  Ak^c?  +  kla^  —  Ak^ac  —  Akik^a^c  <  0  (4.23) 

This  equation  then  factors  to  the  form 

{k^a  +  k^a^  —  2fcgc)^  <  0  (4.24) 


which  is  a  contradiction  to 


{k^a  +  k^o?  —  2kecY  >  0.  (4-25) 

Thus,  the  system  has  no  oscillatory  behavior. 
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The  value  for  A  comes  from  equation  (4.22),  by  way  of  the  quadratic  formula 
A±  =  jf  ^  fQj-  values  of  a  or  the  kinetic  constants,  then  A+  >  0 

and  the  steady  state  curve  has  divergent  behavior.  Now  4ac  =  4:{k4kjyiy4  +  k^kjy^  + 
2/C4/C72/22/4  +  ki^k^y^)  from  equation  (4.22).  Substituting  the  steady  state  points  in 
for  1/1,  2/2,  and  1/4  given  as  equations  (4.15),  (4.16),  and  (4.18)  and  simplifying  gives 
the  value  4ac  =  Now  if  k^c^  -  k^ca  +  2kWc  -  k^a^  <  0, 

then  there  exists  divergent  behavior.  By  factoring,  ^50(0  —  a)  +  k4a^{2c  —  a)  <  0, 
but  from  equation  (4.17),  c  >  a, for  all  values  of  c  and  a,  so  the  steady  state  curve 
is  stable  for  all  c  and  a. 

The  zero  eigenvalue,  presented  in  equation  (4.21),  introdiices  the  possibility 
that  the  system  represented  by  equations  (4.4)- (4.6)  may  be  rediiced  to  a  system 
of  two  equations.  Reducing  the  system  can  provide  further  understanding  of  the 
models  steady  state  behavior. 

Note  from  equations  (4.4)-(4.6),  that  ^  ~  ^  This  implies  that 

^  +  ^  —  ^=0,  or  more  appropriately  that  2/4  —  2/2  +  yi  =  C2  where 

C2  =  2/4(0)  -  2/2(0)  -f  2/i(0).  (4.26) 

By  substituting  2/4  =  C2  +  2/2  —  2/i  5  equations  (4.4)  and  (4.6)  become 


k4yiy2  +  k^ic  -  2/2)  (4.27) 

k4yiy2  +  (^5  +  ke){c  -  2/2)  -  fc72/2(c2  +  y2  -  yi)  (4.28) 

At  steady  state,  values  can  be  located  by  setting  ^  =  0  and  ^  =  0.  The  resulting 
equations  are 


dyi 

(it 

dy2 

(it 
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0 


(4.29) 

(4.30) 


=  -Kym  +  h{.c  -  y2) 


0  =  -k^y^y2  +  {h  h){c  -  y2)  -  k^y2{c2  +  y2  -  yi) 


Solve  equations  (4.29)  and  (4.30)  for  yiand  ?/4  to  get  the  steady  state  concentrations 


2/1 

2/2 


-^( 

1 

2k7ki 


2ckik4  +  k^ke  +  ^702^4  +  k'jk^  — 
k^kf,  +  k'fC2k4  +  k'jk^  — 

•  {-k^ke  -  k7C2k4  -  k^k^  +  y/r^) 


(4.31) 

(4.32) 


where  r2  =  kfk^  +  2k^kQk^C2  +  2k4kQk'fk^  +  k^c^kl  +  2k^C2k4k^  +  k'^k'l  +  Ak^k^k^c  + 
Ak^k^k^c.  Equations  (4.31)  and  (4.32)  involve  no  free  parameters,  so  by  reducing 
the  system  of  equations  (4.4)-(4.6)  to  two  equations  given  by  (4.27)  and  (4.28),  the 
steady  state  line  becomes  a  steady  state  point.  Plugging  in  the  initial  conditions 
2/1(0)  =  50,  2/2(0)  =  100,  2/3(0)  =  0  ,  2/4(0)  =  0  so  that  c  =  2/2(0)  +  2/3(0)  =  100 
and  C2  =  2/4(0)  —  2/2(0)  +  2/1(0)  =  —50  where  k^  =  2,  k^  =  1,  k^  =  2,  and  /C7  =  1, 

gives  r2  =  13  025,  such  that  2/1  =  0.456  36,  2/2  =  52.  282,  2/3  =  47.  718,  and  2/4  =  1- 
825  6.  These  results  are  consistent  with  the  results  produced  previously  by  the  three 
equation  system. 

The  Jacobian  of  the  system  (4.27)-(4.28)  takes  the  form 


J 


-k4y2  -k^yi  -  k^ 

-k4y2  +  k’jy2  -k^yi  -k^-k^-  ^702  -  2k^'y2  +  /c7yi 


(4.33) 


Then  solving  for  the  characteristic  equation,  det|AJ  —  J|  =  0 


det 


A  +  k4y2  k^yi  +  k^ 

km  -  km  A  +  km  +  k5  +  kG  +  km  +  2km  -  kjyi 


0  (4.34) 
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gives  the  polynomial  +  6iA  +  bo  =  0.  where  bi  =  k^yi  +  +  ke  +  /C7C2  +  2/C72/2  — 

kjUi  +  /C42/2,  and  bo  =  k4^key2  +  ^702^4^/2  +  ‘^kjk^y^  +  k'jk^y2-  The  eigenvalues  are 
given  by  A  =  \{—bi  ±  Now,  A  is  very  complicated  and  although 

an  analytic  evaluation  may  be  possible,  simplifying  A  is  difficult.  One  reason  for 
this  is  due  to  the  negative  value  within  equation  (4.26).  The  negative  values  allow 
the  possibility  for  C2  <  0,  so  C2  cannot  be  used  to  simplify  A.  For  this  reason,  a 
numerical  evaluation  is  more  practical.  The  eigenvalues  corresponding  to  the  steady 
state  point  evaluated  above  are  A+  =  —55.967  and  A_  =  —106.61.  The  negative 
eigenvalues  correspond  to  stable  behavior  within  a  neighborhood  of  the  steady  state 
point  (2/i,y2,y3,l/4)ss  =  (0.456  36,52.282,47.718,1.8256). 

This  system  observed  steady  state  behavior  in  the  form  of  a  stable  steady 
state  curve  when  evaluated  as  a  system  of  three  differential  equations  and  as  a 
stable  steady  state  point  when  the  system  was  reduced  to  two  equations.  The  three 
equation  system  had  a  single  zero  eigenvalue  corresponding  to  a  free  parameter  in 
the  steady  state.  The  zero  eigenvalue  and  free  parameter  were  eliminated  when  the 
system  was  reduced  and  although  the  original  system  consisting  of  four  equations  was 
not  evaluated,  it  is  likely  that  it  has  two  zero  eigenvalues  corresponding  to  two  free 
parameters.  It  is  worth  noting  that  the  behavior  of  the  system  was  actually  easier 
to  analytically  evaluate  when  reduced  to  three  equations,  than  when  reduced  to  two 
equations.  Although  the  system  was  analyzed  to  look  for  different  forms  of  behavior 
near  the  steady  state  curve,  specifically  oscillatory  convergence  or  divergence,  the 
analysis  shows  that  only  asymptotically  stable  behavior  exists  within  this  system  for 
positive  real  values  of  kinetic  constants  and  the  parameter  a. 

This  small  Glucose  model  shows  that  some  systems  may  exhibit  stable  behavior 
over  a  curve,  yet  perturbing  concentrations  within  the  system  lead  them  to  flow  to 
a  different  concentration  on  that  curve.  A  system  may  stay  within  a  region  of 
the  steady  state  curve  in  the  following  way.  With  constant  small  perturbations 
of  one  type,  the  steady  state  concentrations  may  flow  down  that  curve.  Another 
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type  of  perturbation  may  cause  the  concentration  to  flow  back  up  the  curve,  thus 
regulating  the  system.  Hypothetically,  when  modeling  a  cell,  the  steady  state  curve 
may  contain  concentrations  that  are  toxic.  A  large  series  of  the  same  type  of 
perturbation  could  lead  the  steady  state  concentrations  to  these  toxic  levels  and  kill 
the  cell.  Therefore,  the  concentrations  introduced  into  this  smaller  system  must  be 
regulated  through  perturbations  that  are  controlled  elsewhere  in  the  larger  model. 

4.6.2  Large  Glucose  Model.  The  smaller  Glucose  model  is  contained  in 
a  larger  system  that  is  listed  in  Appendix  A  and  described  in  Chapter  3,  Section 
3.7.  This  larger  model  incorporates  the  smaller  system  as  both  the  second  and 
eighth  of  ten  total  reactions.  In  this  model,  reactions  one,  three,  and  ten  are  one 
way  reactions  that  degrade  the  reactant  concentrations  of  all  previous  reactions  to 
a  steady  state  value  of  zero.  Thus,  the  second  reaction,  denoting  the  small  glucose 
model,  shows  no  interesting  results  when  analyzed  as  a  part  of  the  larger  model. 
That  is,  if  reaction  three  is  included,  then  the  steady  state  of  reaction  two  is  zero 
for  all  of  its  substrates.  However,  since  reaction  eight  is  the  small  model  and  all 
reactions  between  reaction  three  and  reaction  ten  are  two  way  reactions,  if  reaction 
ten  is  omitted  from  the  system,  the  large  system  can  be  run  to  steady  state  and 
the  small  Glucose  model  can  be  analyzed  as  a  part  of  the  larger  model.  In  this 
case,  reaction  eight  acts  as  the  smaller  model  and  must  cooperate  with  the  two  way 
reactions,  four  through  nine. 

Analysis  of  how  the  smaller  system  relates  to  the  larger  system  takes  place  in 
steps.  By  making  some  assumptions  about  the  larger  glucose  model,  reaction  eight 
acts  just  the  same  as  the  small  glucose  model  already  evaluated.  The  assumptions 
that  give  the  larger  glucose  model  the  characteristics  of  the  smaller  glucose  model 
are  as  follows.  First,  by  setting  the  concentrations  of  all  the  substrates  (includ¬ 
ing  catalyst  transition  reactants)  except  3-phosphoglycerate  and  2-phosphoglycerate 
equal  to  zero,  system  changes  will  be  initiated  only  from  concentrations  of  reactants 
involved  in  reaction  eight.  Then,  by  setting  the  concentrations  of  the  enzymes. 
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Figure  4.16.  The  large  glucose  model  with  assumptions  that  make  it  behave  as 
though  it  were  the  small  glucose  model. 

Phosphoglycerate  Kinase  of  reaction  seven  and  Enolase  of  reaction  nine,  equal  to 
zero,  the  reacting  species  are  restricted  to  only  reaction  eight.  This  can  also  be  ac¬ 
complished  by  setting  the  concentration  of  the  coenzyme  Mg++  equal  to  zero.  The 
resulting  model  acts  just  the  same  as  the  small  glucose  model.  Setting  the  initial 
concentration  of  3-Phosphoglycerate  =  50,  Phosphoglycerate  Mutase  =  100,  the  cat¬ 
alyst  transition  =  0.00001,  and  2-Phosphoglycerate  =  0.00001,  with  forward  reaction 
kinetic  constants  k27  and  k29  equal  to  2  and  reverse  reaction  kinetic  constants  k28 
and  k30  equal  to  1,  the  systems  behavior  reaches  steady  state  quickly.  The  Figure 
4.16  indicates  how  the  substrates  and  enzyme  change  concentrations  through  time 
to  reach  steady  state.  Much  of  the  concentration  of  these  reactants  is  eaten  up  by 
the  catalyst  transition,  which  at  steady  state  equals  47.59. 

By  allowing  the  initial  concentration  of  a-D-Glucose  to  be  a  nonzero  value, 
and  by  setting  the  initial  concentrations  of  the  enzymes,  Phosphoglycerate  Kinase  of 
reaction  seven  and  Enolase  of  reaction  nine,  equal  to  nonzero  values,  reaction  eight 
now  interacts  with  the  larger  system.  Evaluating  this  system  with  small  incremental 
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increases  in  a-D-Glucose,  Phosphoglycerate  Kinase,  and  Enolase  individually,  allows 
for  comparison  between  the  large  glucose  model  and  the  small  glucose  model. 

To  begin,  the  model  was  given  incremental  increases  in  the  concentration  of 
a-D-Glucose.  With  the  concentration  of  Phosphoglycerate  Kinase  =  1  and  the 
concentration  of  Enolase  =  0.  Figure  4.17  describes  how  concentrations  within  the 
small  system  move  to  steady  state  when  a-D-Glucose  =  20  initially.  In  this  case, 
reaction  eight  quickly  reaches  steady  state  just  as  it  did  in  the  previous  system. 
However,  as  this  is  occurring,  concentrations  of  reactants  originating  from  a-D- 
Glucose  are  flowing  through  the  system.  At  a  later  time,  these  concentrations  react 
to  create  3-Phosphoglycerate  and  2-Phosphoglycerate.  As  greater  amounts  of  these 
substrates  are  created,  a  larger  amount  of  the  catalyst  transition  in  reaction  eight  is 
also  created,  leading  to  a  depleting  of  the  amount  of  the  enzyme  Phosphoglycerate 
Mutase.  A  small  decrease  in  3-Phosphoglycerate  also  occurs  rapidly  because  the 
small  amount  of  Phosphoglycerate  Kinase  introduced  into  the  system,  coupled  with  a 
small  reverse  reaction,  tries  to  balance  the  large  concentration  of  3-Phosphoglycerate 
with  the  small  concentration  of  the  substrate,  1,3-Biphosphoglycerate. 

The  model  was  then  given  incremental  increases  in  Phosphoglycerate  Kinase, 
with  the  concentrations  of  a-D-Glucose  =  0  and  Enolase  =  0.  Figure  4.18  shows 
that  as  the  concentration  of  Phosphoglycerate  Kinase  increases,  a  reverse  reaction 
quickly  takes  effect  in  the  model  by  removing  concentrations  of  3-Phosphoglycerate 
to  balance  its  concentration  with  the  concentration  of  1,3-Bisphosphoglycerate. 


Finally,  the  model  was  given  incremental  increases  in  Enolase,  with  the  con¬ 
centrations  of  a-D-Glucose  =  0  and  Phosphoglycerate  Kinase  =  0.  Figure  4.19 
shows  that  as  the  concentration  of  Enolase  increases,  a  forward  reaction  takes  effect 
in  the  model  and  removes  concentrations  of  2-Phosphoglycerate  as  it  is  created  until 
a  balance  in  concentration  between  2-Phosphoglycerate  and  Phosphoenolpyruvate 
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when  alpha-D-Glucose  =  20 


Steady  state 


Figure  4.17.  Incremental  increase  in  the  concentration  of  a- D- Glucose 


when  Phosphoglycerate  Kinase  =  20 


Steady  state 


Figure  4.18.  Incremental  increase  in  the  concentration  of  Phosphoglycerate  Kinase 
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when  Enolase  =  5 


Figure  4.19.  Incremental  increase  in  the  concentration  of  Enolase 

is  created.  Thus,  more  3-Phosphogly cerate  and  Phosphoglycerate  Mutase  must  be 
used  up  before  measurable  values  of  2-Phosphoglycerate  begin  to  appear. 

Analysis  of  the  larger  system,  omitting  reaction  ten,  shows  similarities  to  the 
smaller  Glucose  model.  In  fact,  for  the  values  checked  though  simulation,  the  steady 
state  curve  for  the  larger  glucose  model  is  virtually  equal  to  the  one  associated  with 
the  small  glucose  model.  In  the  smaller  system,  if  a  small  perturbation  is  created, 
the  model  will  not  return  to  the  original  steady  state  point,  but  will  go  to  a  new 
steady  state  point  on  the  steady  state  line.  In  the  larger  model,  a  small  perturbation 
that  either  increases  or  decreases  one  or  many  reactant  concentrations  to  disrupt  the 
system,  also  causes  each  reactant  to  move  to  a  different  steady  state  value  defined  on 
a  steady  state  curve.  In  this  case,  because  the  system  is  closed,  the  small  increase  or 
decrease  in  concentration  must  be  preserved.  So  the  system  resets  its  steady  state 
to  account  for  the  change  in  concentrations.  This  change  in  steady  state  can  be 
tracked  with  respect  to  each  participating  reactant.  The  resulting  changes  of  steady 
state  act  as  a  stable  steady  state  cmve. 
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4-7  Gepasi 

The  Gepasi  program,  outlined  in  Chapters  2  and  3,  was  used  to  verify  the 
results  created  in  Matlab.  The  benefit  of  using  a  simiilator  such  as  Gepasi  is  that  it 
allows  for  a  very  friendly  user  interface,  so  verifying  results  is  both  quick  and  easy. 
In  this  section,  the  small  and  large  glucose  models  created  and  executed  in  Matlab, 
will  be  compared  to  simulations  of  these  models  executed  in  Gepasi.  By  evaluating 
a  system  using  both  Matlab  and  Gepasi,  the  results  can  be  compared  to  check  for 
errors  in  the  program. 

The  Matlab  results  at  5000  time  steps  for  the  small  ghicose  model  with  ki¬ 
netic  constant  values  of  k4=2,  k5=l,  k6=2,  and  k7=l  and  initial  conditions  of  3- 
Phosphoglycerate  =  50,  Phosphoglycerate  Mutase  =  100,  Catalyst  transition  =  0, 
and  2- Phosphoglycerate  =  0.00001  are  given  by  the  following. 


3-Phosphoglycerate 
Phosphoglycerate  Mutase 
Catalyst  transition 
2-Phosphoglycerate 


0.4563, 

(4.35) 

52.2817, 

(4.36) 

47.  718, 

(4.37) 

1.8254 

(4.38) 

Substituting  the  same  conditions  into  the  Gepasi  simulator  gives  the  same 
values  at  5000  time  steps. 

[  Yl]  =  4.563564e-001  mM, 

[  Y2]  =  5.228177e+001  mM, 

[  Y3]  =  4.771824e+001  mM, 

[  Y4]  =  1.825425e+000  mM, 

Although  this  technique  isn’t  generally  needed  for  small  systems,  it  is  valuable 
for  bigger  systems.  Consider  the  output  after  5  time  steps  for  the  large  glucose 
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model  ^1  listed  in  Appendix  A.  The  quick  key  gives  a  general  description  of  the 
concentrations  the  output  is  describing.  For  a  more  detailed  understanding,  see 
Appendix  A. 

1 _ QUICK  KEY _ 

7o 

7  CO-ENZYME  =  y(l),  Mg++ 

7  ENZYMES  =  y(2)  -  y(ll) 

7  CATALYST  TRANSITION  =  y(12)  -  y(21) 

7  ENERGY  =  y(22)  -  Y(26) 

7  SUBSTRATES  =  y(27)  -  y(37) 

7 _ 

ans  = 

Columns  1  through  7 

68.4342  999.9999  999.9977  999.9828  999.9132  998.7863  979.5891 
Columns  8  through  14 

80.0914  995.6312  994.6608  993.6994  0.0002  0.0023  0.0172 
Columns  15  through  21 

0.0868  1.2137  20.4109  919.9086  4.3688  5.3393  6.3006 
Columns  22  through  28 

100.4887  973.2848  1.4305  1.4305  978.1586  0.0000  0.0000 
Columns  29  through  35 

0.0000  0.0000  0.0035  0.0005  0.0001  0.0050  0.0001 
Columns  36  through  37 
0.0000  42.2363 

In  this  large  glucose  model,  at  5  seconds,  the  Matlab  results  are  that  ATP 
is  y(22)  =  100.4887  and  ADP  is  y(23)  =  973.2848.  This  relates  quite  well  to  the 
output  from  the  same  model  written  in  Gepasi.  The  output  from  Gepasi,  given  next 
indicates  that  at  5  seconds,  ATP  =  100.19  and  ADP  =  973.5.  The  small  amount  of 
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error  comes  in  part  from  the  different  numerical  analysis  techniques  used  in  the  two 
systems. 


RESULTS  OF  INTEGRATION  (after  5.00e+000  s) 

[  alphaDg]  initial  =  5 . OOOOOOe+002  mM,  final  =  7.673559e-012  mM 
[  ATP]  initial  =  1 . OOOOOOe+003  mM,  final  =  1 .001906e+002  mM 
[  ENZl]  initial  =  1 . OOOOOOe+003  mM,  final  =  6 . 832746e+001  mM 
[  ENZ2]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 . 999999e+002  mM 
[  XI]  initial  =  1 . OOOOOOe-005  mM,  final  =  1 . 575943e-004  mM 
[  alphaDgGp]  initial  =  1 . OOOOOOe-005  mM,  final  =  9 .488750e-007  mM 
[  ADP]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 .735061e+002  mM 
[  ENZ3]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 . 999976e+002  mM 
[  X2]  initial  =  1 . OOOOOOe-005  mM,  final  =  2 . 371148e-003  mM 
[  BataDf6p]  initial  =  1 . OOOOOOe-005  mM,  final  =  1 . 208948e-009  mM 
[  ENZ4]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 . 999823e+002  mM 
[  X3]  initial  =  1 . OOOOOOe-005  mM,  final  =  1 . 771586e-002  mM 
[  BataDflSb]  initial  =  1 . OOOOOOe-005  mM,  final  =  4.742581e-005  mM 
[  ENZ5]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 . 999110e+002  mM 
[  X4]  initial  =  1 . OOOOOOe-005  mM,  final  =  8. 900037e-002  mM 
[  G3p]  initial  =  1 . OOOOOOe-005  mM,  final  =  3 .455546e-003  mM 
[  Dihp]  initial  =  1 . OOOOOOe-005  mM,  final  =  4. 933821e-004  mM 
[  ENZ6]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 . 987871e+002  mM 
[  X5]  initial  =  1 .OOOOOOe-005  mM,  final  =  1 . 212868e+000  mM 
[  NAD]  initial  =  1 . OOOOOOe+003  mM,  final  =  1 .435234e+000  mM 
[  P]  initial  =  1 . OOOOOOe+003  mM,  final  =  1 .435234e+000  mM 
[  ENZ7]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 . 795344e+002  mM 
[  X6]  initial  =  1 .OOOOOOe-005  mM,  final  =  2.046566e+001  mM 
[  Bisl3]  initial  =  1 .OOOOOOe-005  mM,  final  =  5 . 796210e-005  mM 
[  NADH]  initial  =  1 . OOOOOOe-005  mM,  final  =  9.780991e+002  mM 
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[  ENZ8]  initial  =  1 . OOOOOOe+003  mM,  final  =  8.004378e+001  mM 
[  X7]  initial  =  1 .OOOOOOe-005  mM,  final  =  9. 199562e+002  mM 
[  Phos3]  initial  =  1 .OOOOOOe-005  mM,  final  =  5 . 017286e-003  mM 
[  ENZ9]  initial  =  1 . OOOOOOe+003  mM,  final  =  9. 956060e+002  mM 
[  X8]  initial  =  1 .OOOOOOe-005  mM,  final  =  4.394014e+000  mM 
[  Phos2]  initial  =  1 .OOOOOOe-005  mM,  final  =  9 . 054872e-005  mM 
[  ENZIO]  initial  =  1 .OOOOOOe+003  mM,  final  =  9 . 946307e+002  mM 
[  X9]  initial  =  1 .OOOOOOe-005  mM,  final  =  5 . 369277e+000  mM 
[  Phos]  initial  =  1 . OOOOOOe-005  mM,  final  =  1 . 131147e-007  mM 
[  ENZll]  initial  =  1 . OOOOOOe+003  mM,  final  =  9 . 936708e+002  mM 
[  XIO]  initial  =  1 . OOOOOOe-005  mM,  final  =  6 .329212e+000  mM 
[  Pyruvate]  initial  =  1 . OOOOOOe-005  mM,  final  =  4 . 204531e+001  mM 

By  performing  an  independent  calculation,  a  numerical  algorithm  can  be  par¬ 
tially  validated.  The  Gepasi  simulator  was  a  great  asset  in  getting  a  better  under¬ 
standing  of  how  the  larger  system  should  behave. 

J^.8  S^immany 

Chapter  4  gave  insight  into  how  small  perturbations  in  model  parameters  and 
variables  affect  the  system  as  a  whole.  Perturbations  on  the  linear  model  showed 
variations  of  steady  state  behavior  within  a  linear  system.  Perturbations  on  the 
Brusselator  showed  problems  associated  with  analyzing  a  linearized  model.  Per¬ 
turbations  of  both  reactant  concentrations  and  constant  parameters  in  the  Schnack- 
enberg  model  showed  great  sensitivity  in  the  form  of  large  amplihcation  that  may 
have  broad  affects  on  other  systems  within  a  larger  model.  Hasty’s  bistable  system 
showed  how  different  magnitudes  of  perturbations  of  reactant  concentrations  can 
permanently  change  the  behavior  of  a  system.  The  small  Glucose  model  extended 
stable  steady  states  points  to  stable  steady  state  curves  and  showed  how  perturba¬ 
tions  in  reactant  concentrations  can  lead  to  different  steady  states  within  the  curve, 
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while  perturbations  within  parameters  can  change  the  shape  of  the  curve.  The 
larger  Glucose  model  then  extended  the  analysis  of  small  systems.  Incorporating 
the  small  glucose  model  into  the  larger  Glucose  model  leads  to  insight  on  how  to 
incorporate  small  systems  into  larger  ones.  Furthermore,  it  shows  that  doing  so  can 
lead  to  a  better  understanding  of  the  larger  system.  Finally,  the  glucose  models 
were  applied  to  the  model  simulator  Gepasi  as  an  example  of  how  this  analysis  can 
be  verified  using  user  friendly  simulators  available  over  the  internet. 
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V.  Summary  and  Conclusion 


5. 1  Summarxj 

The  goal  of  this  thesis  was  to  understand  the  underlying  behavior  of  small 
models  and  their  relation  to  larger  systems  in  an  effort  to  better  understand  how  the 
Air  Force  should  construct  intracellular  models  to  assist  in  ftitiire  toxicology  studies. 
The  approach  taken  in  this  thesis  was  to  create  a  system  of  reaction  equations, 
and  then  transform  that  system  into  a  system  of  differential  eqiiations  that  could 
be  solved.  Several  models  were  derived  using  this  technique  and  then  analyzed. 
New  models  were  created  including  two  models  of  transcription  and  translation,  a 
large  model  of  glycolysis,  and  a  small  component  model  of  the  large  ghicose  model. 
Then,  a  number  of  the  previously  derived  models  as  well  as  additional  models  were 
evaluated  using  a  sensitivity  analysis  in  an  effort  to  identify  how  their  behavior  may 
change  with  small  perturbations  in  both  initial  conditions  and  constant  parameters. 
This  analysis  technique  was  applied  to  the  newly  created  small  glycolysis  model  and 
incorporated  into  the  larger  model  in  an  effort  to  understand  how  the  behavior  of 
smaller  models  affects  larger  systems. 

5.2  Conclusion 

Reaching  the  goal  of  this  thesis  has  led  to  a  number  of  conclusions.  First,  the 
technique  of  modeling  a  physical  system  using  biochemical  reaction  equations  and 
then  converting  the  reaction  equations  to  a  system  of  ordinary  differential  equations 
is  appropriate  for  smaller  systems.  Second,  analysis  of  models  using  the  technique  of 
linearizing  about  steady  state  and  then  determining  the  stability  of  the  linear  model 
shows  that  it  is  capable  of  displaying  interesting  behavior  in  small  systems.  Third, 
this  technique  can  be  used  to  analyze  larger  systems. 

The  technique  of  using  reaction  equations  to  model  transcription,  translation, 
and  metabolism  was  a  result  of  reviewing  the  modeling  techniques  presented  by 
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Hasty  et  al.  [14],  [13],  Schnackenberg  [30],  and  Prigogine  et  al.[27].  By  deriving 
results  concerning  existing  models  that  involved  reaction  eqiiations,  it  was  hoped 
that  more  understanding  would  be  gained  as  to  how  to  develop  new  models. 

Analyzing  a  model  involved  starting  with  the  reaction  eqiiations  and  trans¬ 
forming  them  into  a  system  of  ordinary  differential  eqiiations  that  were  then  used 
to  determine  the  steady  state  behavior  of  the  system.  Examples  of  steady  state 
behavior  in  the  form  of  both  stable  behavior  and  unstable  behavior  were  evaluated 
by  deriving  linear  systems  from  nonlinear  systems  using  linearization.  Hypothet¬ 
ical  systems  including  the  Brusselator  equation  and  Schnackenbergs  eqiiation  were 
linearized  and  then  the  behavior  near  the  steady  state  was  determined  from  the  eigen¬ 
values  of  the  linearized  system.  A  small  system  with  multiple  steady  state  points 
was  evaluated  to  determine  how  to  find  the  steady  states  of  the  system.  This  knowl¬ 
edge  was  incorporated  in  solving  for  the  steady  state  points  of  two  real-life  cellular 
systems  developed  by  Hasty  et  al.  Using  the  information  gathered  from  evaluat¬ 
ing  these  systems,  two  small  reaction  models  involving  transcription  and  translation 
were  presented  and  derived.  A  third  reaction  model  involving  the  glycolysis  pathway 
was  created  in  the  form  of  a  system  of  35  ordinary  differential  equations  and  then 
solved  using  a  Matlab  script. 

To  discover  how  concentrations  of  reactants  behave  under  uncertain  conditions, 
some  of  the  derived  models  were  analyzed  graphically  to  trace  the  behavior  of  the 
systems.  The  linear  system  was  evaluated  first.  Results  showed  behavior  as  it 
relates  to  stable  and  unstable  steady  state  points.  This  analysis  was  expanded  to 
nonlinear  systems  by  evaluating  the  Brusselator  system.  Under  different  conditions, 
the  Busselator  showed  both  convergence  to  a  limit  cycle  and  oscillatory  convergence 
to  a  steady  state  point.  This  analysis  was  then  used  to  test  the  variability  of  various 
models  with  respect  to  changing  initial  conditions  and  uncertain  constant  parame¬ 
ters.  By  analyzing  the  Schnackenberg  model,  under  changing  initial  conditions  with 
one  set  of  parameters,  it  was  shown  that  the  convergence  to  a  limit  cycle  might  be 
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more  dramatic  when  the  initial  conditions  fall  outside  the  limit  cycle.  Changing 
the  constant  parameters  slightly  showed  great  variability  in  how  the  model  behaved 
as  well.  A  small  change  in  constant  parameters  gave  rise  to  oscillatory  stability 
that  showed  great  sensitivity  to  initial  conditions.  A  slightly  larger  change  in  con¬ 
stant  parameters  created  straight  line  stability.  Next,  by  considering  Hasty  et  al.’s 
model,  small  changes  in  initial  conditions  led  concentrations  to  move  to  different 
stable  steady  state  points.  A  steady  state  analysis  was  performed  on  Hasty  et 
al.’s  model  in  the  form  of  a  bifurcation  diagram  to  determine  steady  state  behavior 
with  respect  to  uncertain  constant  parameters.  With  the  knowledge  gained  from 
evaluating  these  models,  a  small  subcomponent  of  the  newly  created  ghicose  model 
was  evaluated.  After  reducing  the  system  from  four  eqiiations  to  three,  evaluation 
showed  it  had  multi-stability  in  the  form  of  a  stable  steady  state  curve  for  every 
variation  of  parameters  and  initial  conditions  within  reasonable  model  boundaries. 
Rirther  reduction  to  a  system  of  two  equations  removed  a  free  parameter  and  allowed 
the  reduced  system  to  be  evaluated  for  a  single  steady  state  point. 

In  an  effort  to  understand  how  small  models  relate  to  larger  models,  the  small 
glucose  model  was  incorporated  into  the  larger  glycolysis  model.  The  small  sys¬ 
tem  was  evaluated  under  various  conditions  of  the  larger  model.  The  evaluation 
examined  changes  that  took  place  in  the  smaller  model  and  how  it  was  affected  by 
conditions  imposed  by  the  larger  system. 

Altogether,  the  information  gathered  within  this  thesis  presents  a  good  basis  for 
understanding  how  to  evaluate  cellular  systems.  It  also  provides  an  introduction  for 
how  the  Air  Force  should  construct  intracellular  models  to  assist  in  future  toxicology 
studies. 


5. 3  Rexommendations 

Further  research  into  the  construction  and  evaluation  of  cellular  models  would 
benefit  Air  Force  toxicology  studies.  It  remains  to  be  determined  how  many  types 
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of  behavior,  such  as  limit  cycle  convergence,  are  attained.  Continued  analysis  of  the 
models  described  in  this  thesis,  as  well  as  other  models,  should  lead  to  understanding 
of  how  these  types  of  behavior  come  about. 

Continued  analysis  can  be  performed  by  evaluating  Hasty  et  al.’s  [14]  model 
without  the  assumption  of  fast  reactions.  Hasty  et  al.’s  model  might  also  be  ex¬ 
tended  to  include  promoter  molecules  and  evaluated  to  see  how  the  behavior  changes. 
Also,  new  models  ^1  and  #2  can  be  evaluated  numerically  to  determine  whether 
or  not  they  display  interesting  behavior.  Further  insight  into  how  to  build  models 
displaying  interesting  behavior  may  be  gained  by  incorporating  models  displaying 
different  types  of  behavior  into  a  larger  model.  For  example,  by  incorporating  a 
sensitive  model,  such  as  the  Schnackenberg  model,  into  a  system  that  involves  miilti- 
stability,  such  as  Hasty  et  al.’s  model,  a  small  perturbation  of  reactant  concentration 
in  the  sensitive  system  may  be  amplified  in  a  way  that  makes  reactant  concentrations 
of  the  multi-stable  model  move  to  another  stable  steady  state.  Also,  by  coupling  the 
model  described  by  Chen  et  al.  [7]  with  the  glucose  model,  transcription,  transla¬ 
tion,  and  metabolism  can  be  tied  together.  Chen’s  model  can  be  used  to  regulate 
the  10  enzymes  of  the  gluoce  model  thereby  controling  the  prodution  of  ATP.  Also, 
techniques  for  solving  for  the  steady  state  of  larger  systems  should  be  explored.  This 
might  be  done  by  comparing  the  speed  of  numerical  techniques  with  simply  running 
the  program  to  steady  state. 

Models  not  included  in  this  thesis  that  should  be  considered  for  future  evalua¬ 
tion  include  stochastic  models  and  time  delay  models.  Finally,  hypothetical  cellular 
models  involving  toxicology  issues  should  be  created,  separated  into  smaller  mod¬ 
els,  analyzed,  and  numerically  evaluated  to  give  insight  into  what  types  of  behavior 
should  be  studied  in  greater  detail. 
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Appendix  A.  The  Large  Glucose  Model 

Two  models  of  Glucose  are  included  here.  The  first  is  a  user  friendly  version 
that  allows  the  user  to  understand  and  track  all  the  reactants  within  the  model. 
Furthermore,  the  output  is  displayed  in  an  easy  to  read  manner  and  all  inputs  into 
the  program  can  be  made  easily.  This  is  an  appropriate  model  to  use  when  viewing 
the  system  as  a  whole.  In  this  case,  it  is  used  to  evaluate  the  changes  in  ATP,  ADP, 
and  the  intermediate  concentrations.  Due  to  Matlabs  difficult  programming  nature, 
evaluating  individual  aspects  of  the  system  requires  a  more  complicated  version  of 
the  program.  This  more  complex  version  is  included  as  a  way  to  view  individiial 
components  of  the  model,  and  make  changing  the  system  an  easier  task.  In  this 
version,  the  last  equation  removed  as  a  way  to  test  the  small  Glucose  model  in  a 
larger  system  under  conditions  of  steady  state. 

For  the  convenience  of  the  user,  each  large  Glucose  model  was  written  as  two 
Matlab  programs.  The  first  program  displays  the  output  and  allows  the  user  to 
change  initial  conditions.  It  then  calls  the  next  program.  The  second  program 
contains  the  large  Glucose  model.  All  the  detail  about  these  programs  is  included 
as  comments  in  the  programs. 

A.l  Large  Glucose  Model  #1 
A.  1.1  Program  1. 

7o  yO  is  the  initial  concentrations  of  the  system. 

y0=[1000;  1000;  1000;  1000;  1000;  1000;  1000;  1000;  1000;  1000;  1000; 


00001; 

.00001; 

.00001; 

.00001;  .00001; 

.00001; 

.00001;  .00001; 

00001; 

.00001; 

1000;  1000;  1000;  1000; 

.00001; 

500;  .00001; 

00001; 

.00001; 

.00001; 

.00001;  .00001; 

.00001; 

.00001;  .00001; 

00001]  ; 
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tsp=[0,5]  7o  this  is  the  time  interval  for  the  calculation 

[t, y]=ode23s(’ Glucose ’ ,tsp,yO) ;  %  Here  the  Glucose  program  is  called 

n=max(size(y)) ; 

y(n,:)  7  prints  out  the  concentration  at  the  final  time  step 
7  Evaluating  ATP  and  ADP  concentrations 

yT  =  y(:,12)  +  y(:,14)  +  y(:,18)  +  y(:,21)  +  y(:,22)  +  y(:,23); 
yp  =  [y(:,12)  y(:,14)  y(:,18)  y(:,21)  y(:,22)  y(:,23)  yT] 
plot(t ,yp) 

7  In  the  plot,  ATP  =  y(22)  and  ADP  =  y(23)  concentration  changes 
7  involve  intermediate  transition  steps  y(12),  y(14) ,  y(18) ,  and 
7  y(21) .  The  total  concentration  of  these  values  though  time  is 
7  equal  to  the  initial  concentrations.  Notice  that  y(12)  uses  up 
7  almost  all  ATP  immediately,  then  gets  almost  all  used  up  by  the 
7  time  y(21)  is  being  created. 

A.  1.2  Program  2. 

function  [dydt]  =  dydt(t,y) 
n=max(size(y)) ; 
dydt=ones (37 , 1) ; 

7  There  are  37  different  molecules  that  are  changing  in  this  system. 
7  Each  reaction  has  a  reaction  rate  coefficient  associated  with  it. 

7  The  reaction  rate  coefficients  are  given  in  this  code  by  kl-k37. 

7  The  forward  reaction  rate  coefficients  are  3  and  the  backwards 
7  reaction  rate  coefficients  are  1.  In  this  model,  all  reaction 
7  rate  coefficients  are  fixed.  In  future  models,  they  may  be 
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7o  considered  variables  or  functions.  NOTE:  If  there  is  no  backward 
7  reaction,  then  the  reaction  is  a  one-way  reaction, 
kl  =  3;  7REACTI0N  #1 
k2  =  1; 

k3  =  3;  7REACTIQN  #2 
k4  =  3;  7REACTIQN  #3 
k5  =  1; 

k6  =  3;  7REACTI0N  #4 
k7  =  1; 

k8  =  3;  7REACTI0N  #5 
k9  =  1; 

klO  =  3;  7REACTI0N  #6 
kll  =  3;  7REACTI0N  #7 
kl2  =  1; 

kl3  =  3;  7REACTI0N  #8 
kl4  =  1; 

kl5  =  3;  7REACTI0N  #9 
kl6  =  1; 

kl7  =  3;  7REACTI0N  #10 
kl8  =  1; 

kl9  =  3;  7REACTI0N  #11 
k20  =  1; 

k21  =  3;  7REACTiaN  #12 
k22  =  1; 

k23  =  3;  7REACTI0N  #13 
k24  =  1; 

k25  =  3;  7REACTI0N  #14 
k26  =  1; 
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k27  =  3;  7oREACTION  #15 
k28  =  1; 

k29  =  3;  7oREACTION  #16 
k30  =  1; 

k31  =  3;  7oREACTIQN  #17 
k32  =  1; 

k33  =  3;  BEACTIQN  #18 
k34  =  1; 

k35  =  3;  '/REACTION  #19 
k36  =  1; 

k37  =  3;  '/REACTION  #20 

'/  Each  reaction  is  dependent  the  molecules  formed  by  previous 

'/  reactions.  Therefore,  all  reaction  in  this  system  are 

'/  interconnected  through  time.  This  interconnection  makes  modeling 

'/  the  system  as  a  whole,  a  very  difficult  process.  However,  first 

'/  looking  at  each  reaction  individually,  then  incorporating 

'/  the  equations  for  the  individual  reactions  into  the  complete  model 

'/  makes  this  model  easier  to  construct. 

'/  There  are  10  individual  reactions  in  the  very  simple  biochemical 
'/  model  of  the  metabolic  pathway.  Glycolysis.  However,  each 
'/  reaction  is  a  catalytic  reaction.  This  being  the  case,  each 
'/  reaction  can  be  split  into  two  reactions  that  more  accurately 
'/  describe  the  catalytic  reaction.  SEE:  Mathematics  Applied  to 
'/  Deterministic  Problems  in  the  Natural  Sciences  by  Lin,  Segel, 

'/  and  Handleman.  This  changes  the  10  reactions  into  20 

'/  slightly  more  complex  reactions  involving  hypothetical  molecular 
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7o  states  of  catalytic  transition.  Each  of  the  20  reactions  has  a 
7  reaction  rate  depending  on  the  concentration  of  molecules  entering 
7  into  the  system  and  the  concentration  of  those  leaving.  Before 
7  looking  at  the  system  as  a  whole,  the  individual  reaction  systems 
7  are  constructed  here . 

7 _ QUICK  KEY _ 

7  CO-ENZYME  =  y(l),  Mg++ 

7  ENZYMES  =  y(2)  -  y(ll) 

7  CATALYST  TRANSITION  =  y(12)  -  y(21) 

7  ENERGY  =  y(22)  -  Y(26) 

7  SUBSTRATES  =  y(27)  -  y(37) 

7 _ 


7  Mg++  +  Hexokinase  + 

7  alpha-D-Glucose  +  ATP  (k2)<=  =>(kl)  Catalyst_transition_l 
7  is  equivalent  to 

eql  =  -kl*y(27)*y(22)*y(l)*y(2)  +  k2*y(12) ; 

7  Mg++  +  Hexokinase  + 

7  Catalyst_transition_l  =>(k3)  alpha-D-Glucose-6-phosphate  +  ADP 
7  is  equivalent  to 
eq2  =  -k3*y(12); 

7  Phosphoglucoisomerase  + 

7  alpha-D-Glucose-6-phosphate  (k5)<=  =>(k4)  Catalyst_transition_2 

7  is  equivalent  to 

eq3  =  -k4*y(28)*y(3)  +  k5*y(13) ; 
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7o  Phosphoglucoisomerase  + 

7  Catalyst_transitioii_2  (k7)<=  =>(k6)  beta-D-Fructose-6-phosphate 

7  is  equivalent  to 

eq4  =  -k6*y(13)  +  k7*y(29)*y(3) ; 

7  Mg++  +  Phosphofructokinase  + 

7  beta-D-Fructose-6-phosphate  +  ATP  (k9)<=  =>(k8) 

Catalyst_transition_3 

7  is  equivalent  to 

eq5  =  -k8*y(29)*y(22)*y(l)*y(4)  +  k9*y(14) ; 

7  Mg++  +  Phosphofructokinase  + 

7  Catalyst_transition_3  =>(klO)  beta-D-Fructose-1 ,6-bisphosphate  +  ADP 
7  is  equivalent  to 
eq6  =  -kl0*y(14); 

7  Aldolase  + 

7  beta-D-Fructose-1 , 6-bisphosphate  (kl2)<=  =>(kll) 

Catalyst_transition_4 

7  is  equivalent  to 

eq7  =  -kll*y(30)*y(5)  +  kl2*y(15); 


7  Aldolase  + 

7  Catalyst_transition_4  (kl4)<=  =>(kl3)  Glyceraldehyde-3-ph.oshate 

+  Dihydroxyacetone  phosphate 

7  is  equivalent  to 

eq8  =  -kl3*y(15)  +  kl4*y(31)*y(32)*y(5) ; 


A-6 


7o  Triose  phosphate  isomerase  + 

7o  Dihydroxy  ace  tone  phosphate  (kl6)<=  =>(kl5)  Catalyst_transition_5 

7  is  equivalent  to 

eq9  =  -kl5*y(32)*y(6)  +  kl6*y(16); 

7  Triose  phosphate  isomerase  + 

7  Catalyst_transition_5  (kl8)<=  =>(kl7)  Glyceraldehyde-3-phoshate 

7  is  equivalent  to 

eqlO  =  -kl7*y(16)  +  kl8*y(31)*y(6) ; 

7  Glyceraldehyde-3-phosphate-dehydrogenase  + 

7  (NAD+  +  Pi)  +  Glyceraldehyde-3-phoshate  (k20)<=  =>(kl9) 

Catalyst_transition_6 

7  is  equivalent  to 

eqll  =  -kl9*y(31)*y(25)*y(24)*y(7)  +  k20*y(17) ; 

7  Glyceraldehyde-3-phosphate-dehydrogenase  + 

7  Catalyst_transition_6  (k22)<=  =>(k21)  NADH  +  1,3-Bisphosphoglycerate 
7  is  equivalent  to 

eql2  =  -k21*y(17)  +  k22*y(33)*y(26)*y(7) ; 

7  Mg++  +  Phosphoglycerate_kinase  + 

7  1,3-Bisphosphoglycerate  +  ADP  (k24)<=  =>(k23)  Catalyst_transition_7 
7  is  equivalent  to 

eql3  =  -k23*y(33)*y(23)*y(l)*y(8)  +  k24*y(18); 

7  Mg++  +  Phosphoglycerate_kinase  + 

7  Catalyst_transition_7  (k26)<=  =>(k25)  3-Phosphoglycerate  +  ATP 
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7o  is  equivalent  to 

eql4  =  -k25*y(18)  +  k26*y(34)*y(22)*y(l)*y (8) ; 

7o  Phosphoglycerate_mutase  + 

7  3-Phosphoglycerate  (k28)  <=  =>  (k27)  Catalyst_transition_8 

7  is  equivalent  to 

eql5  =  -k27*y(34)*y(9)  +  k28*y(19); 


7  Phosphoglycerate_mutase  + 

7  Catalyst_transition_7  (k30)<=  =>(k29)  2-Phosphoglycerate 
7  is  equivalent  to 
eql6  =  -k29*y(19)  +  k30*y(35)*y(9) ; 

7  Mg++  +  Enolase  + 

7  2-Phosphoglycerate  (k32)<=  =>(k31)  Catalyst_transition_9 
7  is  equivalent  to 

eql7  =  -k31*y(35)*y(l)*y(10)  +  k32*y(20); 

7  Mg++  +  Enolase  + 

7  Catalyst_transition_9  (k34)<=  =>(k33)  Phophoenolpyruvate 
7  is  equivalent  to 

eql8  =  -k33*y(20)  +  k34*y(36)*y(l)*y(10) ; 

7  Mg++  +  Pyruvate_Kinase  + 

7  Phophoenolpyruvate  +  ADP  (k36)<=  =>(k35)  Catalyst_transition_10 
7  is  equivalent  to 

eql9  =  -k35*y(36)*y(23)*y(l)*y(ll)  +  k36*y(21) ; 
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7o  Mg++  +  Pyruvate_Kiiiase  + 

7  Catalyst_transitioii_10  =>(k37)  Pyruvate  +  ATP 
7  is  equivalent  to 
eq20  =  -k37*y(21); 

7  Now  each  of  the  37  molecules  change  with  respect  to  time. 

7  Each  of  the  equations  is  the  reaction  rate 
7  written  with  respect  to  a  single  molecule  on  the  left  hand 
7  side  of  the  biochemical  equation.  The  negative  of  that 
7  equation  will  give  the  reaction  rate  written  with  respect 
7  to  a  single  molecule  on  the  right  hand  side  of  the  biochemical 
7  equation.  To  find  the  reaction  rate  of  a  specific  molecule 
7  over  the  entire  system,  add  each  equation  where  the  molecule 
7  appears  on  the  left  hand  side  of  the  biochemical  equation. 

7  Of  course,  if  the  molecule  appears  on  the  right  hand  side 
7  of  the  biochemical  equation,  add  the  negative  of  the  above 
7  equations. 

dydt(l)  =  eql  -  eq2  +  eq5  -  eq6  +  eql3  -  eql4 
+  eql7  -  eql8  +  eql9  -  eq20; 

7  Change  of  Mg++  with  respect  to  time 
dydt(2)  =  eql  -  eq2; 

7  Change  of  Hexokinase  with  respect  to  time 
dydt(3)  =  eq3  -  eq4; 

7  Change  of  Phosphoglucoisomerase  with  respect  to  time 
dydt(4)  =  eq5  -  eq6; 

7  Change  of  Phosphofructokinase  with  respect  to  time 
dydt(5)  =  eq7  -  eq8; 
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7o  Change  of  Aldolase  with  respect  to  time 
dydt(6)  =  eq9  -  eqlO; 

7o  Change  of  Triose  phosphate  isomerase  with  respect  to  time 
dydt(7)  =  eqll  -  eql2; 

7  Change  of  Glyceraldehyde-3-phosphate-dehydrogenase 
7  with  respect  to  time 
dydt(8)  =  eql3  -  eql4; 

7  Change  of  Phosphoglycerate  kinase  with  respect  to  time 
dydt(9)  =  eql5  -  eql6; 

7  Change  of  Phosphoglycerate  mutase  with  respect  to  time 
dydt(lO)  =  eql7  -  eql8; 

7  Change  of  Enolase  with  respect  to  time 
dydt(ll)  =  eql9  -  eq20; 

7  Change  of  Pyruvate  kinase  with  respect  to  time 
dydt(12)  =  -eql  +  eq2; 

7  Change  of  Catalyst_transition_l  with  respect  to  time 
dydt(13)  =  -eq3  +  eq4; 

7  Change  of  Catalyst_transition_2  with  respect  to  time 
dydt(14)  =  -eq5  +  eq6; 

7  Change  of  Catalyst_transition_3  with  respect  to  time 
dydt(15)  =  -eq7  +  eq8; 

7  Change  of  Catalyst_transition_4  with  respect  to  time 
dydt(16)  =  -eq9  +  eqlO; 

7  Change  of  Catalyst_transition_5  with  respect  to  time 
dydt(17)  =  -eqll  +  eql2; 

7  Change  of  Catalyst_transition_6  with  respect  to  time 
dydt(18)  =  -eql3  +  eql4; 

7  Change  of  Catalyst_transition_7  with  respect  to  time 
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dydt (19) 
7o  Change 
dydt (20) 
7o  Change 
dydt (21) 
7  Change 
dydt (22) 
7  Change 
dydt (23) 
7  Change 
dydt (24) 
7  Change 
dydt (25) 
7  Change 
dydt (26) 
7  Change 
dydt (27) 
7  Change 
dydt (28) 
7  Change 
dydt (29) 
7  Change 
dydt (30) 
7  Change 
dydt (31) 
7  Change 
dydt (32) 
7  Change 


=  -eql5  +  eql6; 

of  Catalyst_transition_8  with  respect  to  time 
=  -eql7  +  eql8; 

of  Catalyst_transition_9  with  respect  to  time 
=  -eql9  +  eq20; 

of  Catalyst_transition_10  with  respect  to  time 
=  eql  +  eq5  -  eql4  -  eq20; 
of  ATP  with  respect  to  time 
=  -eq2  -  eq6  +  eql3  +  eql9; 
of  ADP  with  respect  to  time 
=  eqll; 

of  NAD+  with  respect  to  time 
=  eqll; 

of  Pi  with  respect  to  time 
=  -eql2; 

of  NADH  with  respect  to  time 
=  eql; 

of  alpha-D-Glucose  with  respect  to  time 
=  -eq2  +  eq3; 

of  alpha-D-Glucose-6-phosphate  with  respect  to  time 
=  -eq4  +  eq5; 

of  beta-D-Fructose-6-phosphate  with  respect  to  time 
=  -eq6  +  eq7 ; 

of  beta-D-Fructose-l,6-bisphosphate  with  respect  to  time 
=  -eq8  -  eqlO  +  eqll; 

of  Glyceraldehyde-3-phoshate  with  respect  to  time 
=  -eq8  +  eq9; 

of  Dihydroxyacetone  phosphate  with  respect  to  time 
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dydt (33) 
7o  Change 
dydt (34) 
7o  Change 
dydt (35) 
7  Change 
dydt (36) 
7  Change 
dydt (37) 
7  Change 


=  -eql2  +  eql3; 

of  1,3-Bisphosphogly cerate  with  respect  to  time 
=  -eql4  +  eql5; 

of  3-Phosphoglycerate  with  respect  to  time 
=  -eql6  +  eql7 ; 

of  2-Phosphoglycerate  with  respect  to  time 
=  -eql8  +  eql9; 

of  Phophoenolpyruvate  with  respect  to  time 
=  -eq20; 

of  Pyruvate  with  respect  to  time 


Figure  A.l  is  a  graph  that  describes  how  ATP  and  ADP  move  through  the 
larger  model.  This  larger  model  is  simply  an  example  of  how  small  systems  can 
be  incorporated  into  larger  systems.  In  this  example,  at  time  t  =  0,  ATP  has  a 
concentration  of  1000  units,  ADP  has  a  concentration  of  1000  units,  and  the  four 
intermediate  steps  have  a  concentration  of  .0001  units. 


It  is  clear  that  throughout  this  system,  for  the  first  five  time  steps,  ATP,  ADP, 
and  the  intermediate  steps  change  concentrations  quite  rapidly.  However,  it  is  also 
clear  by  the  line  that  represents  the  total  concentrations  of  all  the  reactants,  that 
no  reactants  diverge.  This  is  a  property  that  held  true  for  the  smaller  system  and 
is  also  true  for  the  larger  system. 

As  time  moves  out  to  5000  time  steps,  it  is  clear  that  the  reactants  still  do  not 
diverge.  Stable  systems  such  as  these  are  what  biological  systems  rely  on  to  survive. 
Through  these  systems,  they  may  control  their  to  and  begin  to  settle  down  into  a 
form  that 


A. 2  Large  Glucose  Model  ^2 
A.  2.1  Program  1. 


A- 12 


Time 


Figure  A.l.  Concentration  graph  of  ATP,  ADP  and  associated  catalyst  transitions 
in  the  large  glucose  model 


Time 


Figure  A.2.  Concentration  graph  of  ATP,  ADP  and  associated  catalyst  transitions 
in  the  large  glucose  model 
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7o  Initial  concentrations 

y0=[500;  10;  1000;  100;  1000;  .00001;  .00001;  1000;  .00001;  .00001 
1000;  .00001;  .00001;  1000;  .00001;  .00001;  .00001;  1000;  .00001; 
1000;  1000;  .00001;  1000;  .00001;  .00001;  10;  .00001;  100.00001; 
1000;  .00001;  .00001;  10;  .00001;  .00001]; 

tsp=  [0,5000]  7  plug  in  the  time  interval 
[t,y]=ode23s(’Glucose_tapered.’  ,tsp,y0) ; 
m=max(size(y)) ; 

7  Evaluating  reaction  concentrations 
c=  y(: ,29)  +  y(: ,30) ; 

yp  =  [y(:,28)  y(:,29)  y(:,30)  y(:,31),  c] 

7  plots  the  concentrations  of  each  reactant  in  reaction  8 
7  with  respect  to  time 

7plot(t ,yp) 

k31  =  1.1;  7REACTI0N  #17 
k32  =  1; 

k33  =  1.1;  7REACTI0N  #18 
k34  =  1; 
cO  =  c(l,l); 

step  =  m;  7  Input  the  time  step  (an  integer  between  1  and  n) 
alpha  =  y (step, 29); 
time_step  =  t(step); 

7 
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7o  Checks  the  steady  state  value 
y_28=(c0*k32  -  alpha*k32)/(alpha>i'k31) ; 
y_31=(c0*k33  -  alpha*k33)/(alpha*k34) ; 

[y(m,:),  t(m)]  7  prints  out  the  concentration  at  the  final  time  step 
steady_state  =  [time_step  y_28  alpha  cO-alpha  y_31] 

7 

7  plots  how  the  concentrations  of  the  reactants  move  through  time. 

7  plot3(y(: ,28),  y(:,29),  y(:,31)) 

A. 2. 2  Program  2. 

function  [dydt]  =  dydt(t,y) 
n=max(size(y) ) ; 
dydt=ones (34 , 1 ) ; 

7  This  is  the  number  of  elements 
7  in  this  pathway  that  go  through 
7  a  change  over  time. 

7  These  are  the  reaction  rate  coefficients 
7  The  reaction  rates  are  taking  place  with 
7  rate  coefficient  concentrations  of  millimoles/sec. 

kl  =  3;  7REACTI0N  #1 
k2  =  .01; 

k3  =  3;  7REACTI0N  #2 
k4  =  3;  7REACTI0N  #3 
k5  =  .01; 

k6  =  3;  7REACTI0N  #4 
k7  =  .01; 

k8  =  3;  7REACTI0N  #5 
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k9  =  .01; 

klO  =  3;  7oREACTION  #6 
kll  =  3;  7oREACTION  #7 
kl2  =  .01; 

kl3  =  3;  7oREACTIQN  #8 
kl4  =  .01; 

kl5  =  3;  “/oREACTIQN  #9 
kl6  =  .01; 

kl7  =  3;  7REACTI0N  #10 
kl8  =  .01; 

kl9  =  3;  °/oREACTION  #11 
k20  =  .01; 

k21  =  3;  y„REACTION  #12 
k22  =  .01; 

k23  =  3;  “/oREACTION  #13 
k24  =  .01; 

k25  =  30;  y„REACTION  #14 
k26  =  .0001; 

k27  =  1.1;  “/oREACTION  #15 
k28  =  1; 

k29  =  3;  y„REACTION  #16 
k30  =  .01; 

k31  =  1.1;  y„REACTIQN  #17 
k32  =  1; 

k33  =  1.1;  yREACTION  #18 
k34  =  1; 

y  k35  =  3;  yREACTION  #19  ycommented  out 
y  k36  =  1 ;  ycoimnented  out 


A- 16 


7o  k37  =  3;  7oREACTION  #20  7conmented  out 


A  =  -kl*y(l)*y(2)*y(4)*y(5)  +  k2*y(6); 

B  =  -k3*y(6); 

C  =  -k4*y(7)*y(8)  +  k5*y(9) ; 

D  =  -k6*y(9)  +  k7*y(8)*y(10) ; 

E  =  -k8*y(10)*y(2)*y(4)*y(ll)  +  k9*y(12) ; 

F  =  -kl0*y(12); 

G  =  -kll*y(13)*y(14)  +  kl2*y(15); 

H  =  -kl3*y(15)  +  kl4*y(14)*y(16)*y(17) ; 

I  =  -kl5*y(17)*y(18)  +  kl6*y(19); 

J  =  -kl7*y(19)  +  kl8*y(18)*y(16); 

K  =  -kl9*y(16)*y(20)*y(21)*y(23)  +  k20*y(24) ; 

L  =  -k21*y(24)  +  k22*y(22)*y(23)*y(25) ; 

M  =  -k23*y(25)*y(3)*y(4)*y(26)  +  k24*y(27) ; 

N  =  -k25*y(27)  +  k26*y (2) *y (4) *y(26) *y (28) ; 

0  =  -k27*y(28)*y(29)  +  k28*y(30); 

P  =  -k29*y(30)  +  k30*y(29)*y(31) ; 

Q  =  -k31*y(31)*y(4)*y(32)  +  k32*y(33); 

R  =  -k33*y(33)  +  k34*y (4) *y (32) *y (34) ; 

7  S  =  -k35*y(34)*y(3)*y(4)*y(35)  +  k36*y(36) ;  7conmiented  out 
7  T  =  -k37*y(36);  7commented  out 

dydt ( 1 )  =  A ; 

dydt(2)  =A+E-N;  7-T;  7values  commented  out 
dydt (3)  =  -B  -F  +  M;  7  +  S;  7values  commented  out 
dydt(4)  =A-B+E-F+M-N+q-R;  7+S-T; 

7  values  of  dydt (4)  are  commented  out 

dydt(5)  =  A  -  B; 
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dydt(6)  =  -A  +  B; 
dydt(7)  =  -B  +  C; 
dydt(8)  =  C  -  D; 
dydtO)  =  -C  +  D; 
dydt(lO)  =  -D  +  E; 
dydt(ll)  =  E  -  F; 
dydt(12)  =  -E  +  F; 

dydt(13)  =  -F  +  G; 

dydt(14)  =  G  -  H; 
dydt(15)  =  -G  +  H; 

dydt(16)  =  -H  -  J  + 

dydt(17)  =  -H  +  I; 

dydtdS)  =  I  -  J; 
dydt(19)  =  -I  +  J; 

dydt(20)  =  K; 

dydt(21)  =  K; 

dydt(22)  =  -L; 

dydt(23)  =  K  -  L; 
dydt(24)  =  -K  +  L; 

dydt(25)  =  -L  +  M; 

dydt(26)  =  M  -  N; 
dydt(27)  =  -M  +  N; 

dydt(28)  =  -N  +  0; 

dydt(29)  =  0  -  P; 
dydt(30)  =  -0  +  P; 

dydt(31)  =  -P  +  Q; 

dydt(32)  =  Q  -  R; 
dydt(33)  =  -Q  +  R; 


dydt(34)  =  -R;  % 
7o  dydt(35)  =  S  - 
7  dydt(36)  =  -S 
7  dydt(37)  =  -T; 


,  +  S;  7values  connnented  out 
T;  7conimeiited  out 
+  T ;  7conimented  out 
7coifflneiited  out 
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